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CHAPTER  1 
INTRODUCTION 


1.1  BACKGROUND 

Projectile  impact  and  penetration  problems  can  include  the  attack  of  a  military 
target  with  penetrating  weapons,  propellant  embedded  anchors,  storm  launched  debris 
impacting  a  structure,  etc.  Analysis  of  these  problems  requires  knowledge  of  (1)  the 
impact  conditions  (velocity  and  projectile  orientation  relative  to  the  target)  (2) 
characteristics  of  the  projectile  (geometry,  mass,  and  strength)  and  (3)  the  properties  of 
the  target  (dimensions  and  mechanical  properties).  Within  the  limits  of  the  conventional 
impact  velocities  of  interest,  less  than  about  1  km/s,  the  depth  of  penetration  increases 
when  the  kinetic  energy  of  the  projectile  increases.  This  can  be  accomplished  by 
increasing  either  the  impact  velocity  or  the  projectile  mass.  Deformation  of  the  projectile 
can  also  influence  penetrability.  A  heavy,  essentially  non-deformable  projectile,  such  as 
an  armor-piercing  (AP)  projectile,  will  penetrate  deeper  into  a  target  than  a  deformable 
projectile.  An  oblique  impact  results  in  a  highly  asymmetrical  stress  distribution  on  the 
projectile  that  can  cause  it  to  either  breakup  or  result  in  a  lower  penetration  depth  than  a 
normal  impact.  Strong  geomaterials  such  as  concrete  or  rock  have  greater  mass, 
modulus,  and  strength  characteristics  than  soil  and  will  therefore  have  a  greater  resistance 
to  penetration.  The  interaction  of  these  three  categories  of  parameters  is  complex  and  at 
times  leads  to  difficulties  in  interpreting  and  analyzing  impact  and  penetration  events. 

In  general,  depending  on  the  complex  interaction  of  all  the  variables,  any  of  the 
following  situations  could  prevail  when  a  projectile  impacts  a  target.  The  projectile  could 
either  (1)  break  up,  deform  significantly,  and  be  “defeated”  by  the  target,  (2)  ricochet,  (3) 
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initially  penetrate  into  the  target  and  then  broach,  (4)  perforate  the  target,  or  (5)  penetrate 
and  come  to  rest  in  the  target  (Figure  1.1).  The  response  of  target  materials  to  projectile 
impact  will  depend  on  many  variables  such  as  the  material  properties,  impact  velocity, 
projectile  shape,  and  target  size.  Material  may  be  ejected  from  the  target  due  to  spalling 
or  scabbing.  Spalling  is  a  tensile  failure  that  results  from  the  reflection  of  the  initial 
compressive  wave  from  the  rear  surface  of  a  finite  thickness  target.  Scabbing  is  similar  in 
appearance  to  spalling  but  is  formed  by  the  fracturing  and  break  up  of  target  material  due 
to  large  deformations.  These  conditions  may  exist  imder  high  impulse  loadings  due  to 
intense  impacts  into  materials  that  are  stronger  in  compression  than  in  tension,  like  most 
geomaterials.  Under  the  impact  and  subsequent  loading  conditions  of  interest  here,  most 
of  the  target  response  near  free  surfaces  is  described  by  scabbing  and  fracture. 

1 .2  ANALYSIS  TECHNIQUES 

Projectile  penetration  events  are  generally  analyzed  by  three  approaches: 
empirical  curve  fitting,  analytic  models,  and  numerical  methods.  The  selection  of  a 
predictive  method  is  usually  dictated  by  the  level  to  which  the  penetration  event  is  to  be 
analyzed.  If  only  depth  of  penetration  is  of  interest,  then  any  of  the  three  predictive 
methods  may  be  used.  If  the  depth  of  penetration  as  well  as  the  state  of  stress  within  the 
target  is  of  interest,  some  analytical  methods  and  any  of  the  numerical  methods  can  be 
used.  Detailed  description  of  the  stresses  and  the  deformations  of  the  target  require  the 
use  of  the  numerical  approaches. 

The  degree  to  which  the  target  materials  are  characterized  generally  increases  with 
each  approach.  In-depth  analyses  must  account  for  the  geometry  of  the  target,  large  fimte 
strains  and  deflections,  strain  rate  effects,  work  hardening,  heating  or  fiictional  effects, 
and  the  initiation  and  propagation  of  fracture  (Jonas  and  Zukas  1978).  The  empirical 
techniques  (Young  1972,  Bernard  1977,  and  Pahl  1989)  involve  curve  fitting  of 
penetration  test  data  to  relate  depth  of  penetration  and  projectile  deceleration  to  projectile 
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Figiire  1.1.  Projectile  response  after  impact. 
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geometry,  impact  velocity,  and  target  material  type.  The  empirical  relations  are  the 
simplest  to  apply  and  are  generally  reliable  within  the  range  of  conditions  covered  by  the 
data  set  on  which  they  are  based.  The  projectile  is  assumed  to  be  rigid  and  the  target 
material  is  characterized  by  simple  engineering  parameters,  e.g.,  density,  unconfined 
compressive  strength,  penetrability  index  (such  as  S-number;  Yormg  1972).  The 
trajectory  of  the  projectile  is  assumed  to  follow  a  straight  line.  The  empirical  techmques 
offer  no  insight  to  the  physics  or  mechanics  of  penetration  and  the  response  of  the  target 
material. 

The  analytic  models  offer  a  somewhat  more  fundamental  approach  in  that  they  are 
based  on  the  conservation  and  balance  laws  of  continuum  mechanics.  Many  of  these 
models  use  dynamic  cavity  expansion  techniques  (CET)  to  determine  the  resistance  of  the 
target  material  to  projectile  penetration  (Bernard  and  Creighton  1976,  Forrestal  1986, 
Forrestal  and  Luk  1992,  and  Forrestal  and  Tzou  1997).  The  models  are  often  restricted 
by  the  assumptions  made  on  the  target  material  motion  (one-dimensional  spherical  or 
cylindrical  geometries)  and  stress  response  (compressive)  in  order  to  simplify  the 
analysis.  The  simplified  governing  equations  can  usually  be  solved  in  closed  fi’om.  The 
constitutive  properties  of  the  target  materials  required  by  a  given  model  are  determined 
from  independent  laboratory  tests.  The  properties  may  include  density,  compressibility, 
strength,  shear  modulus,  etc.  Recently,  CET  has  been  used  to  develop  algorithms  to 
predict  penetration  into  soil  (Forrestal  and  Luk  1992),  rock  (Forrestal  1986),  and  concrete 
(Forrestal  and  Tzou  1997).  A  limitation  on  the  CET  is  that  it  is  only  valid  for  normal 
impact,  penetration  beyond  the  cratering  phase  and  into  the  tunneling  phase,  and 
penetration  into  thick  targets  where  the  back  surface  does  not  mfluence  the  penetration 
process.  The  model  can  be  compared  to  test  data  and  its  predictive  capability  assessed. 
Disagreements  between  the  predicted  and  actual  penetration  event  may  suggest  that  the 
model  is  not  adequate,  but  may  also  suggest  physical  effects  that  can  not  be  accormted  for 
in  the  simplified  analytical  procedure. 
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Another  analytic  technique  is  the  so-called  differential  area  force  law  (DAFL) 
approach.  The  DAFL  type  of  formulation  (Henderson  and  Stephens  1972)  divides  the 
projectile  longitudinally  and  then  circumferentially  into  a  mesh  of  elements,  or 
differential  areas,  on  the  projectile  surface.  A  normal  stress  is  calculated  on  each  element 
using  stress  formulations  that  may  be  empirical  or  analytic  based.  The  stress  is  divided 
by  the  area  of  the  element  and  the  total  force  is  obtained  by  summing  over  all  the 
elements. 

The  most  comprehensive  approach  to  projectile  penetration  problems  is  the 
numerical  approach  using  finite-element  or  finite-difference  wave  propagation  codes. 

The  numerical  methods  solve  the  continuity,  momentum  and  energy  balance  equations  of 
continuum  mechanics  in  conjimction  with  an  appropriate  constitutive  representation  for 
the  target  materials  of  interest.  These  first-principles  techniques  can  use  a  wide  variety  of 
initial  and  boundary  conditions  to  simulate  the  penetration  event.  The  complexity  of  the 
constitutive  material  models  that  can  be  used  with  most  of  the  numerical  methods  has 
little  restriction.  However,  for  problems  concerning  geomaterials,  fairly  simple  material 
flow  (associated  flow  rules  of  plasticity)  and  fracture  models  (reduction  of  the  intact 
properties  of  the  material)  are  usually  employed  due  to  a  lack  of  understanding  of  how 
these  complex  processes  should  be  modeled.  The  material  property  data  required  for  use 
in  the  numerical  methods  must  be  obtained  from  the  appropriate  independent  laboratory 
tests  on  the  target  materials. 

The  material  models  used  in  the  analysis  of  penetration  problems  must 
incorporate  the  physical  phenomena  controlling  the  process.  Ideally,  the  model  should 
account  for  compaction,  cracking,  shear  dilation,  water  migration,  phase  transformation, 
thermal  effects,  inhomogeneity,  etc.  Geomechanical  models  (such  as  the  Prandtl-Reuss 
model,  CAP-type  models,  etc.;  Chen  1982)  that  simulate  elasto-plastic  deformations, 
compaction,  and  failure  are  available,  but  the  implementation  must  be  applicable  to  the 
high-pressure,  sub-millisecond  loadings  that  occur  during  high-velocity  projectile  impact. 
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This  may  require  that  the  model  be  rate-dependent  and  include  an  equation-of-state  to 
define  the  pressure-volume  relation  for  the  very  high,  impulsive  loading.  Fracture  and 
damage  models  based  on  the  micromechanical  approach  that  simulate  the  opening  of 
existing  cracks  and  the  formation  of  new  cracks  have  been  developed,  but  they  have  not 
been  applied  to  penetration  problems  involving  geomaterials  for  several  reasons.  Most  of 
these  models  are  developed  to  simulate  the  response  of  materials  at  low  pressures  under 
static  or  quasi-static  loading  conditions.  Interpretation  of  the  model  parameters  can  be 
obscure  and  difficult  to  obtain  from  laboratory  experiments.  If  implemented  into  a  finite- 
element  code,  computation  time  can  be  prohibitive  due  to  the  large  number  of  cracks  and 
the  volume  of  damaged  material  that  are  produced  during  projectile  penetration.  A  need, 
therefore,  exists  to  improve  the  geologic  material  models  used  to  analyse  penetration  into 
geomaterials  to  accoimt  for  pertinent  phenomena  such  as  brittle  failure,  post-fi'acture 
material  response,  and  pressure  and  rate  effects,  and  to  incorporate  the  models  into 
numerical  analysis  codes. 

As  pointed  out  by  Desai  and  Siriwardane  (1994),  formulation  of  a  viable 
constitutive  law  involves  the  following  steps:  (1)  develop  the  mathematical  formulation, 
(2)  identify  the  significant  parameters,  (3)  determine  the  parameters  fi'om  laboratory 
experiments,  (4)  successful  prediction  of  a  majority  of  observed  laboratory  data  firom 
which  the  model  parameters  were  determined  and  prediction  of  material  response  to  other 
loading  paths,  and  (5)  satisfactory  comparisons  between  numerical  simulations  of 
relevant  bormdary  value  problems  using  the  constitutive  law  and  results  fi'om 
experiments.  The  relevant  boundary  value  problems  of  interest  in  this  research  involve 
the  penetration  and  perforation  of  brittle  geomaterials  by  high-velocity  projectiles. 

Perhaps  a  precursor  to  these  steps  should  be  the  conduct  of  an  experiment,  for  example  a 
beam  test,  or  a  projectile  penetration  experiment,  or  a  blast  loading  in  a  rock-type 
material,  etc.  Results  firom  the  experiments  should  be  observed  carefully  as  to  how  the 
materials  responded.  Does  the  material  have  a  brittle  or  ductile  response,  or  a 
combination  of  the  two?  Does  the  material  contain  many  cracks  and,  if  so,  their  size  and 
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patterns?  Is  the  material  loose  and  friable  after  the  test?  Where  within  the  material  do 
these  conditions  exist?  Step  one  above  should  then  try  to  mimic  these  observations. 

1.3  OBJECTIVES  AND  SCOPE 

The  objectives  of  the  research  reported  herein  were  (1)  develop  a  constitutive 
model  for  brittle  geomaterials  for  numerical  simulation  of  high-velocity  projectile 
penetration  and  perforation  problems,  (2)  determine  the  numerical  values  for  the 
parameters  of  the  model  from  a  series  of  laboratory  material  property  tests,  (3)  implement 
the  model  into  an  existing  large-deformation  finite-element  computer  program,  (4) 
conduct  a  series  of  controlled  laboratory  ballistic  experiments  involving  high-velocity 
penetration  and  perforation  of  plain  concrete  targets,  and  (5)  perform  an  in-depth 
numerical  analysis  of  the  ballistic  experiments  and  evaluate  the  capability  of  the 
constitutive  model  to  capture  the  salient  features  of  the  penetration  and  perforation  of 
brittle  geomaterials. 

Chapter  2  presents  a  review  of  relevant  projectile  impact  and  penetration 
phenomena  with  emphasis  placed  on  brittle  geomaterials.  Chapter  3  contains  the  results 
of  material  property  test  programs  and  includes  relevant  information  on  testing,  analysis 
of  data,  and  stress  paths  of  interest  to  penetration  problems.  Chapter  4  details  the 
development  of  a  constitutive  model  for  brittle  geomaterials  and  the  determination  of  the 
numerical  values  of  the  parameters  of  the  model  for  plain  concretes  used  in  ballistic 
programs.  Implementation  of  the  constitutive  model  into  a  large-deformation  finite- 
element  computer  program  is  presented  in  Chapter  5.  The  ballistic  program  and  the 
comparisons  of  the  results  for  the  penetration  and  perforation  experiments  with 
corresponding  numerical  simulations  are  documented  in  Chapter  6.  Chapter  7  contains 
conclusions  and  recommendations  for  additional  research  related  to  the  topic  of  interest. 
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CHAPTER  2 

PROJECTILE  IMPACT  AND  PENETRATION 

2.1  INTRODUCTION 

The  impact  between  a  projectile  and  a  target  has  been  of  interest  for  hundreds  of 
years.  Historically  the  interest  has  been  towards  military  applications  such  as  the 
penetration  of  hardened  structures  and  armor.  Civilian  interests  in  impact  have  begun  to 
receive  more  study.  These  interests  include  demolition,  transportation  safety, 
crashworthiness  of  vehicles,  weather  borne  impacts,  aircraft  impacts,  body  armor,  erosion 
and  fracture  of  solids,  etc.  (Zukas,  et.al.  1982).  The  phenomena  associated  with  impact 
are  many  and  include  elastic  and  plastic  deformation,  wave  propagation,  fracture  and 
damage,  fiiction,  and,  at  very  high  velocities,  hydrodynamic  flow.  Attempts  have  been 
made  to  classify  impact  by  various  regimes  using  parameters  such  as  striking  velocity 
(Table  2.1;  Zukas,  et.al.  1982)  and  the  strength  of  the  projectile  and  target  and  the  impact 
pressure  (Figure  2. 1 ;  Wilbeck  1985).  According  to  Zukas,  the  problems  of  interest  here 
will  fall  into  the  low  to  intermediate  range  (<1  km/s).  These  impacts  will  result  in 
permanent  damage  where  the  strength  and  compressibility  of  the  materials  are  important. 
Material  loading  and  response  times  are  on  the  order  of  milliseconds  at  low  velocities  and 
microseconds  at  the  higher  velocities.  Based  on  Wilbeck’s  criteria,  the  impacts  will 
generally  involve  conditions  where  the  ratio  of  impact  pressure  to  projectile  strength  is 
less  than  or  approximately  equal  to  one  and  the  ratio  of  impact  pressure  to  target  strength 
is  much  greater  than  one.  The  projectiles  should  sustain  only  slight  deformation,  but  the 
target  will  sustain  significant  deformations. 
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Figure  2.1.  Impact  regimes  based  on  impact  pressure  and  material  strength  (Wilbeck  1985). 
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The  dissipation  of  the  initial  kinetic  energy  of  the  projectile  depends  on  the 
characteristics  of  the  projectile  and  target  and  the  impact  velocity.  In  the  case  of  a  brittle 
target,  much  of  the  projectile’s  energy  will  be  used  to  fragment  and  pulverize  the  target 
(Goldsmith  1960).  The  impact  and  penetration  of  a  projectile  into  a  brittle  geomaterial 
can  be  separated  into  three  possible  phases.  First  is  the  impact  phase  where  the  projectile 
penetrates  the  target  material  only  enough  to  form  an  impact  crater.  The  depth  of  the 
impact  crater  is  typically  only  a  couple  of  projectile  diameters  (calibers)  or  less.  The 
second  phase  is  the  deep  penetration  where  the  projectile  penetrates  beyond  the  impact 
crater  and  begins  to  tunnel  into  the  target  material.  The  tunneling  is  characterized  by  the 
opening  of  a  cylindrical  cavity  by  the  projectile.  Penetration  generally  results  in  the 
embedment  of  the  projectile  in  the  target.  The  third  and  most  complicated  phase  is  the 
perforation  which  is  the  complete  piercing  of  a  target  with  finite  thickness  by  the 
projectile.  This  event  includes  the  formation  of  the  impact  crater,  may  include  the 
tunneling  phase,  and  then  an  exit  condition  that  will  include  the  formation  of  an  exit 
crater.  All  of  these  phases  can  result  in  ejection  of  material  from  the  impact  face  and  the 
back  face  of  the  target  (Bangash  1993)  and  involve  formation  of  cracks,  plastic 
deformation  as  well  as  fragmentation  and  pulverization.  The  three  phases  are  discussed 
in  detail  below. 


2.2  IMPACT 

Impact  can  include  the  collision  of  a  rigid  projectile  with  an  equally  rigid  target,  a 
deforming  projectile  impacting  a  rigid  target,  impact  of  an  essentially  rigid  projectile  with 
a  deforming  target,  or  impact  of  a  deforming  projectile  with  a  deforming  target.  In  some 
cases  the  impact  may  be  followed  by  a  deeper  penetration.  A  rigid  impact  can  be 
analyzed  using  conservation  of  momentum  (Sears,  et.al.,  1976) 


^pVpi  +  =  ^pVpf  +  ^t^tf 


2.1 
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and  conservation  of  energy 


2.2 


where  /Wp  and  w,  are  the  mass  of  the  projectile  and  target,  V^,,-  and  F^,  are  the  initial 
velocity  of  the  projectile  and  target,  and  and  are  the  final  velocity  of  the  projectile 


and  target.  Solving  for  the  final  velocities  yields 
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If  the  target  is  not  only  rigid  with  respect  to  deformation  but  also  with  respect  to 
movement,  F„  =  F^  =  0.0  and  the  final  velocity  of  the  projectile  is  equal  in  magnitude  to 
the  impact  velocity,  but  opposite  in  direction.  This  is  the  simplest  impact  problem  to 
analyze.  It  can  occur  in  several  ways  such  as  very  low  velocity  impacts  between  a 
projectile  and  target  of  strong  materials  and  moderate  velocity  impacts  between  a 
projectile  and  target  of  very  strong  materials.  Once  deformation  begins  to  occur  during 
the  impact,  energy  is  transformed  and  the  analysis  becomes  increasingly  more 
complicated. 

Methods  to  account  for  the  effect  of  deformation  on  the  impact  phase  can  range 
from  the  simple,  such  as  including  Newton’s  coefficient  of  restitution  to  account  for  a 
non-catastrophic  deformation,  to  the  complicated  involving  the  use  of  numerical  methods. 
Two  distinctly  different  modes  of  response  during  impact  of  geomaterials  are  shown 
schematically  in  Figure  2.2,  High-velocity  unpact  into  brittle  geomaterials  will  generally 
result  in  significant  target  deformation  involving  fi-acture,  fi-agmentation  and 
pulverization.  For  non-brittle  materials  the  target  response  is  dominated  by  compaction 
and  shear  flow  and  is  similar  to  that  for  driving  a  pile  into  soil  or  the  bearing  failure  of  a 
shallow  foundation.  Figure  2.3  shows  several  images  taken  during  the  high-velocity 
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Figure  2.2.  Projectile  impact  into  brittle  and  non-brittle  targets. 
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Figure  2.3.  High-velocity  impact  of  a  steel  projectile  into  a  concrete  target  showing  formation  of  the  ejecta  cloud. 
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impact  of  a  steel  projectile  into  a  concrete  target.  The  figure  shows  the  formation  of  the 
ejecta  cloud,  containing  both  the  pulverized  and  fragmented  materials,  as  the  projectile 
enters  the  target.  The  damage  to  a  granite  target  caused  by  the  very  high-velocity  impact 
by  a  steel  projectile  is  shown  in  Figure  2.4  (Ahrens  and  Rubin  1993).  The  figure 
illustrates  and  classifies  the  internal  fractures  radiating  from  the  impact  area  and  the 
highly  fractured  region  near  the  impact  crater.  Concentric,  radial,  spall  and  near  surface 
fractures  are  also  illustrated.  Pressures  near  the  projectile  tip  are  very  high  during  the 
crater  formation  and  cause  the  materials  to  be  pulverized.  The  crater  formation  away 
from  the  projectile  tip  is  primarily  due  to  spall  fracture  and  fragmentation  which  are 
controlled  by  the  shearing  strength  and  tensile  capacity  of  the  target  material.  As  the 
bonds  are  broken  due  to  pulverization  and  fi^gmentation,  the  material  in  these  regions 
begins  to  respond  more  like  a  graniilar  material  than  a  cemented  brittle  material. 

2.3  DEEP  PENETRATION 

Deep  penetration  into  a  semi-infinite  target  follows  the  impact  crater  formation 
and  involves  the  opening  of  a  cylindrical  tunnel  by  the  projectile.  This  phase  is 
dominated  by  severe  compaction  and  high-pressure  shear  flow.  Friction  can  have  some 
effect  during  this  phase,  particularly  as  the  projectile  slows  down.  The  opening  of  the 
tuimel  is  illustrated  in  Figure  2.5.  The  resistance  to  penetration  is  primarily  along  the 
projectile  nose.  As  penetration  progresses,  the  material  surrounding  the  nose  of  the 
projectile  undergoes  severe  compaction  and  shear  flow.  Stresses  near  the  nose  of  the 
projectile  are  extremely  high,  but  dissipate  rapidly  by  80  percent  or  more  only  a  couple  of 
projectile  diameters  away  from  the  penetration  hole. 

The  different  regions  of  material  response  during  deep  penetration  are  described 
in  Figure  2.6  from  the  work  by  Forrestal  (1986)  for  cavity  expansion  analysis.  The  cavity 
expansion  theories  put  forth  by  Forrestal  require  that  the  projectile  penetrate  beyond  the 
impact  phase  and  into  the  tunneling  phase.  The  target  response  during  the  tuimeling 
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Figiire  2.5.  Penetration  of  a  high-velocity  projectile 
into  a  brittle  geomaterial. 
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phase  is  described  as  including  two  to  three  response  regions  (Forrestal  1986).  At  very 
high  velocities,  the  target  response  is  plastic  near  the  penetration  cavity  where  the  stresses 
in  the  material  have  exceeded  its  shear  strength  and  large  deformations  are  occurring. 
Beyond  the  plastic  region,  the  stresses  are  lower  and  the  material  response  is  elastic.  At 
moderate  velocities,  a  cracked  region  is  added  between  the  plastic  and  elastic  region 
where  the  tensile  strength  of  the  material  has  been  exceeded.  At  low  velocities,  where 
penetration  into  the  tuimeling  phase  is  still  possible,  the  plastic  region  is  eliminated  and 
the  material  response  includes  a  cracked  region  followed  by  an  elastic  region. 

Livingston  and  Smith  (1951)  conducted  penetration  experiments  by  air-dropping 
projectiles  into  granite  and  sandstone  targets  from  aircraft  at  different  altitudes.  They 
proposed  that  the  rock  failure  at  impact  is  either  by  plastic  flow  or  by  crushing  and 
fragmentation.  The  failure  of  the  target  material  occurs  in  three  zones  shown  in  Figure 
2.7.  The  zone  of  crushing  occurs  near  the  nose  of  the  projectile  and  is  characterized  by 
the  amount  of  fine  material  found  near  the  nose.  The  amount  of  crushing  depends  on  the 
kinetic  energy  of  the  projectile,  the  sharpness  and  cross-sectional  area  of  the  projectile, 
and  the  properties  of  the  rock.  The  size  of  the  zone  is  greatest  near  the  surface  where  the 
kinetic  energy  of  the  projectile  is  greatest  and  the  confinement  of  the  target  material  is 
lowest.  As  penetration  progresses,  the  size  of  the  zone  decreases  since  the  kinetic  energy 
of  the  projectile  is  decreasing  and  the  confinement  of  the  target  material  is  increasing. 

The  increase  in  confinement  results  in  higher  target  strength.  The  zone  of  shearing  begins 
at  the  zone  of  crushing  and  extends  to  the  zone  of  tension  slabbing.  The  shearing  occurs 
because  of  the  high  compressive  stresses  imparted  to  the  target  material  by  the  projectile, 
but  the  rapid  release  of  these  stresses  as  the  projectile  passes  can  cause  the  material  to 
“burst”  into  the  cavity  behind  the  projectile.  The  zone  of  tension  slabbing  begins  at  the 
zone  of  shearing  and  extends  to  the  surface.  The  extent  of  this  zone  decreases  as 
penetration  progresses.  The  observations  by  Livingston  and  Smith  (1951)  illustrate  the 
importance  of  the  compressibility,  shear  strength,  and  tensile  strength  of  the  target 
material  and  that  the  results  from  a  single  material  response  test,  for  example  an 
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unconfined  compression  test,  cannot  be  used  alone  to  describe  a  target  materials  response 
during  the  penetration  process. 

The  primary  effect  of  friction  between  the  projectile  and  target  in  high-velocity 
projectile  impact  is  near  the  end  of  the  penetration  event.  The  penetration  tunnel  has  a 
tendency  to  close  and  grab  the  cylindrical  aftbody  of  the  projectile  at  this  point  and  bring 
the  projectile  to  a  sudden  stop.  The  influence  of  friction  on  the  penetration  to  this  point  is 
believed  to  be  insignificant  when  compared  to  the  resistance  offered  by  the  target  on  the 
nose  of  the  projectile  as  it  creates  a  tunnel  for  the  projectile  to  go  through. 

2.4  PERFORATION 

Target  perforation  involves  the  formation  of  impact  and  exit  craters,  and  may 
include  a  tunneling  phase  depending  on  the  thickness  of  the  target.  Perforation  will  occur 
in  finite  thickness  targets  where  the  projectile  maintains  sufficient  velocity  throu^  the 
impact  phase  and  the  tunnel  phase,  if  applicable,  to  exit  the  target.  Formation  of  the  exit 
crater  is  a  function  of  the  shear  resistance  and  tensile  strength  of  the  target  material. 
Various  failure  modes  during  perforation  of  finite  thickness  plates  are  shown  in  Figure 
2.8  (Zukas,  et.  al.  1982).  The  perforation  process  may  include  one  of  these  modes  as  the 
dominant  mode,  but  frequently  includes  a  combination  of  several.  In  brittle  geomaterials, 
only  petaling  is  unlikely  to  occur.  Brittle  fracture  will  likely  occur  throughout  the 
perforation  event.  As  the  projectile  nears  an  exit  smface,  radial  cracking,  and  fracture 
may  occur.  The  plug  formation  is  unlikely  to  occur  as  shown  in  Figure  2.8,  but  a  cone  of 
material  may  be  formed  and  pushed  out  from  the  exit  surface  by  the  protruding  projectile. 
Fragmentation  is  likely  if  the  projectile  has  sufficient  velocity  to  pass  through  and  break 
up  the  cone. 

Figure  2.9  shows  the  profiles  of  the  impact  and  exit  craters  for  several 
experiments  using  similar  projectiles  and  impact  velocities,  but  progressively  thinner 


FRAGMENTATION  PETALING- 


Figure  2.8.  Failure  modes  in  impacted  plates  (Backman  1976) 
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Figure  2.9.  Comparison  of  impact  and  exit  craters  for  three  thicknesses  of  unreinforced  concrete  perforated  by  similar  projectiles 
at  similar  impact  velocities. 
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unreinforced  concrete  targets.  The  impact  craters  are  conical  in  shape  and  very  similar  in 
depth  and  width  for  all  experiments.  The  exit  craters  are  also  conical  in  shape,  but  their 
depth  and  width  decreases  as  the  target  thickness  decreases.  Figure  2.10  shows  captured 
images  from  a  high-speed  movie  of  the  back  of  a  concrete  slab  being  perforated  by  a  steel 
projectile.  Radial  cracking  can  be  seen  at  time  tj  as  the  cone  that  was  formed  ahead  of  the 
projectile  is  being  pushed  out  through  the  back  of  the  concrete  slab.  The  projectile  then 
passes  through  the  cone  fracturing  it  into  pieces  of  various  sizes.  Figures  2.9  and  2.10 
indicate  that  the  modes  of  failure  include  combinations  of  brittle  fracture,  radial  cracking, 
and  fragmentation  as  described  in  Figure  2.8. 

Figure  2. 1 1  illustrates  a  hypothesis  for  perforation  of  brittle  materials  and  the 
formation  of  the  exit  crater.  The  projectile  first  creates  the  impact  crater  and  continues  to 
tunnel  through  the  target.  The  tunneling  phase  ends  as  soon  as  the  concentrated  force  on 
the  tip  of  the  protruding  projectile  exceeds  the  resistance  offered  by  a  conical  section  of 
the  target  material  (referred  to  as  the  shear  cone).  The  resistance  offered  by  the  cone  is 
essentially  the  integral  of  the  shearing  strength  of  the  target  material  over  the  surface  area 
of  the  cone.  At  this  point  the  shear  cone  is  formed  and  it  separates  from  the  back  face  of 
the  target  forming  the  exit  crater.  The  pushing  out  of  the  cone  from  the  back  of  the  target 
results  in  a  significant  drop  in  target  resistance  (Pahl  1989).  The  existence  of  a  tunnel 
and  depth  of  the  exit  crater  is  determined  by  the  target  thickness  and  the  strength  of  the 
target  material.  For  very  thick  targets  the  impact  crater  is  followed  by  a  tuimel  and  an 
exit  crater.  As  the  target  becomes  thinner,  the  extent  of  the  tunnel  diminishes.  At  some 
point,  the  tunnel  disappears  completely  and  the  impact  and  exit  meet. 

2.5  SUMMARY  AND  CONCLUSIONS 

A  synopsis  of  the  three  types  of  impact/penetration  events  has  been  presented. 

The  simplest  is  the  rigid  impact  which  can  be  studied  using  momentum  and  energy 
principles.  The  impact  of  a  hard  projectile  into  a  softer  target  will  likely  result  in  the 


25 


d. 

Figure  2.10.  High-speed  movie  images  showing  backface  of  an  unreinforced  concrete  target  being  perforated  by  a  steel  projectile. 
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formation  of  an  impact  crater.  For  brittle  geomaterials,  the  depth  of  the  crater  is 
approximately  two  to  four  projectile  diameters.  The  formation  is  primarily  controlled  by 
the  shear  and  tensile  strength  of  the  target  material.  If  the  impact  conditions  and  the 
extent  of  the  target  are  sufficient  the  formation  of  a  tunnel  phase  may  result.  The 
compressibility  and  strength  of  the  target  determine  the  formation  of  the  tunnel.  For 
finite  thickness  targets,  an  exit  crater  may  form  as  the  projectile  perforates  the  target. 

Much  more  research  has  been  conducted  to  investigate  deep  penetration  than  to 
investigate  impact  and  perforation  of  brittle  geomaterials.  This  illustrates  the  complexity 
of  the  problem  that  requires  modeling  of  the  target  materials  response  to  complex  loading 
paths  in  compression,  tension,  and  shear. 


CHAPTER  3 

MATERIAL  PROPERTIES 


3.1  INTRODUCTION 

In  the  previous  chapter,  the  different  phases  of  impact,  penetration  and  perforation 
were  discussed.  The  material  property  required  for  a  particular  method  of  analysis  is 
dictated  by  the  type  and  complexity  of  the  model  being  used  to  describe  the  target 
material  during  a  penetration  event.  For  the  impact  velocity  range  of  interest  (<  1  km/s), 
the  projectile  is  assumed  to  be  rigid  so  that  only  the  target  response  is  closely  modeled. 
The  empirical,  analytical  and  numerical  methods  of  analysis  each  require  different  levels 
of  target  material  description.  The  empirical  methods  require  a  minimal  amout  of 
information  about  the  target.  The  information  required  by  the  analytic  methods  ranges 
from  a  minimal  amoxmt  similar  to  that  required  by  the  empirical  methods  to  a  more 
detailed  description  requiring  several  specialized  mechanical  property  tests.  The 
numerical  methods  require  detailed  descriptions  of  the  material  response  that  are  typically 
defined  from  results  of  many  specialized  mechanical  property  tests. 

The  detailed  descriptions  of  the  material  properties  are  manifested  in  the 
numerical  methods  through  mathematical  constitutive  models  often  referred  to  as  material 
models.  The  material  models  will  typically  include  a  volumetric  relation,  shear  or 
deviatoric  relation,  and  failure  criteria,  and  may  include  sophisticated  hardening  laws, 
fracture  criteria,  temperature  effects,  and  strain  rate  effects.  The  principal  obstacles  to  the 
use  of  numerical  methods  are  time,  cost,  adequate  description  of  the  problem,  and 
availability  of  adequate  constitutive  models.  Time  and  cost  are  becoming  manageable 
with  the  continued  improvement  of  computational  hardware  and  software.  The  finite- 
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element  codes  use  either  Eulerian  or  Lagrangian  descriptions  to  discretize  the  problem  to 
be  simulated.  Eulerian  codes,  such  as  CTH  (McGlarm,  Thompson,  and  Elrick  1990),  are 
considered  to  be  better  suited  to  problems  involving  large  deformations,  but  sophisticated 
algorithms  are  required  to  propagate  the  various  materials  through  the  fixed  mesh,  and 
simulations  often  require  large  amounts  of  computer  time.  Lagrangian  codes,  such  as 
EPIC  (Johnson,  et  al  1995),  are  more  straightforward,  generally  require  fewer 
computations  per  time  step,  and  the  material  boundaries  are  easily  defined.  Impact 
problems  simulated  using  Lagrangian  codes  require  the  use  of  sliding  interface  algorithms 
to  ensure  separation  of  the  colliding  materials.  As  the  mesh  distortions  become 
significant,  the  accuracy  of  the  computation  decreases  and  the  computational  time  grows 
considerably.  Some  codes  allow  for  the  highly  distorted  regions  to  be  rezoned  either 
manually  or  automatically.  Some  codes  contain  algorithms  that  “erode”  the  highly 
distorted  elements  based  on  the  assumption  that  the  highly  distorted  material  is  no  longer 
contributing  to  the  mechanics  of  the  problem.  Once  the  highly  distorted  elements  have 
been  eroded  the  computation  can  continue  at  a  more  reasonable  pace.  Although  the 
response  of  the  element  material  is  removed  firom  the  calculation,  its  mass  is  conserved  at 
the  nodes.  Material  models  often  found  in  numerical  codes  used  to  simulate  projectile 
penetration  problems  decouple  the  effects  of  volumetric  and  shear  response.  An  equation 
of  state  or  pressure- volume  relation  is  used  to  describe  the  volumetric  response  and 
elastic-plastic  models  are  used  to  describe  the  shear  response. 

The  pressure  levels  developed  during  high-velocity  impact  are  much  higher  than 
those  encountered  in  conventional  loadings.  To  illustrate  the  stress  levels  that  may  be 
encoimtered  during  the  penetration  event,  numerical  simulations  were  made  for  a  steel 
projectile  penetrating  into  a  conventional-strength  (unconfined  compressive  strength  of 
approximately  36  MPa)  concrete  (CSPC)  from  Forrestal,  et.  al.  (1994)  using  the  EPIC 
finite-element  code  (Johnson,  et  al  1994).  The  simulations  were  made  using  an  elastic- 
perfectly  plastic  crushable  solids  material  model.  Figure  3.1  shows  the  description  of  the 
pressure- volume  relation  and  the  strength  relation  for  the  model.  Pressure  is  determined 


in  Figure  3. 1  .a  based  on  the  current  volumetric  strain  and  whether  the  material  is  being 
loaded,  unloaded  or  reloaded.  The  equivalent  stress  o  is  based  on  the  Von  Mises  yield 
criterion  and  is  defined  as 
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o  = 


[(o,  -  *  (o,  -  Oe)^  +  +  ■'?e  +  4)]  3.1 


where  a^,  o^,  Oq,  t„,  T^^e,  and  t^q  are  the  six  stress  components.  A  constant  shear  modulus 
G  is  also  used.  Parameters  used  in  the  model  to  describe  the  CSPC  concrete  are 
summarized  in  Table  3. 1 .  A  comparison  between  the  recommended  material  properties 
for  the  concrete  (Cargile  1998)  and  the  model  fit  is  shown  in  Figure  3.2.  Recommended 
properties  are  based  on  the  results  fi'om  many  triaxial  shear,  hydrostatic  compression,  and 
specialized  (such  as  uniaxial  strain)  experiments  conducted  on  specimens  of  the  concrete. 
They  are  an  interpretation  of  the  experiment  results  to  provide  a  consistant  set  of 
responses  that  can  be  fit  to  mathematical  models  used  in  numerical  simulations.  The 
model  matches  the  hydrostatic  (or  isotropic)  compression  response  of  the  material  fairly 
well,  but,  since  the  ultimate  strength  relation  for  the  model  is  linear,  the  fits  to  the  triaxial 
stress-strain  relations  are  only  approximate. 


Plots  of  Von  Mises  Stress  (VMS)  versus  pressure  (referred  to  as  the  stress  path) 
are  presented  in  Figures  3.3, 3.4  and  3.5  for  several  elements  at  different  ranges  and  at 
depths  below  the  target  surface  of  about  46, 122,  and  635  mm,  respectively,  for  an  impact 
velocity  of  610  m/s.  The  ranges  fi'om  the  impact  point  for  each  figure  are  2.5, 10, 23,  and 
99  mm.  Each  stress  path  shows  an  initial  increase  in  VMS  with  relatively  little  increase 
in  pressure.  The  stress  path  then  follows  the  failure  surface.  The  peak  pressures  at  these 
output  stations  range  fiom  about  840  MPa  to  about  3  MPa.  Several  of  the  stress  paths 
indicate  a  tendency  to  enter  the  tensile  (negative  pressxjre)  region  as  illustrated  by  the 
stress  path  “bumping  into”  the  VMS  axis.  The  calculations  give  some  indication  of  the 
pressure  levels  at  which  material  property  tests  must  be  conducted  to  adequately  capture 
the  material  response  during  a  numerical  simulation  of  the  penetration  event. 


Table  3.1.  Values  for  the  crushable  solids  model  parameters  fit  to  CSPC  concrete. 


mm 

27.0  MPa 

Pcrush 

0.0015 

K, 

10,504.3  MPa 

K2 

-216,518.9  MPa 

K3 

2,604,741  MPa 

Plock 

0.0796 

^lock 

76,000  MPa 

Cl 

37.1  MPa 

C4 

0.49 

s 

1,057.5  MPa 

G 

19,600  MPa 

Density 

2.352  Mg/m^ 

500  I - 1  '  ■  1  . . , .  .  I.  .j  -  ^  ,  — ■  500 


33 


Ddn  *669J)S  lOtUJON  U09|A( 


OdM  *03U9J9|jto  Sfi9j)s  lodpuud 


VON  Kisee  Stress,  osi  (xio^ 


34 


VON  Mtses  Stness.  psi  fxio^  von  Mitee  Stress,  psi  (vio^ 


35 


a.  Range  =  0.1  in.  (2.5  nun) 


b.  Range  =  0.4  in.  (10  mm) 


c.  Range  =  0.9  in.  (23  nun) 


d.  Range  =  3.9  in.  (99  nun) 


Figure  3.4.  Stress-paths  at  a  depth  of  4.8  in.  (122  mm). 
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The  true  properties  of  a  material  are  directly  linked  to  its  structure  at  the  macro¬ 
level,  meso-level,  and  micro-level.  The  properties  are  the  response  of  the  material  to 
mechanical,  physical  or  chemical  influences.  To  determine  the  properties,  an  experiment 
must  be  conducted  on  a  specimen  of  given  size  and  geometry,  under  an  appropriate 
loading  condition  with  specific  boundary  conditions.  The  loading  condition  should  be 
representative  of  the  conditions  that  are  expected  to  exist  in  the  system  as  a  whole  that  is 
to  be  analyzed,  such  as  a  structure,  foundation,  dam,  etc.  (Hordijk,  et.al.  1989). 

The  shear  response  of  geomaterials  is  pressure  dependent.  As  pressure  increases, 
the  ductility  and  strength  of  the  material  changes.  Under  certain  loading  conditions 
volumetric  strains  occur  during  shear  thus  implying  a  coupling  between  the  volumetric 
and  deviatoric  response.  In  the  following  sections,  the  response  of  brittle  geomaterials  to 
hydrostatic  and  deviatoric  states  of  stress  will  be  discussed.  The  influences  of  pressure, 
loading  path,  and  strain  rate  will  be  included. 

3.2  HYDROSTATIC  PRESSURE 

The  pressure- volumetric  strain  response  of  a  geomaterial  is  typically  determined 
by  applying  hydrostatic,  or  equal  all-aroxmd,  pressure  to  a  specimen  of  material  and 
determining  the  resulting  volumetric  strain.  The  pressure  applied  to  the  specimen  is  the 
true  (or  Cauchy)  stress.  The  hydrostatic  pressure  and  volumetric  strain  provide  the  bulk 
response  of  the  material  from  which  the  bulk  modulus  K  is  determined,  usually  as  a 
tangent  to  the  pressure- volumetric  strain  curve.  In  fully-saturated  materials,  the 
volumetric  strain  is  determined  under  drained  conditions  by  monitoring  the  amount  of 
pore  fluid  being  forced  out  of  the  specimen.  Under  imdrained  conditions  the  change  in 
volumetric  strain  is  zero  if  the  compressibility  of  the  pore  fluid  is  neglected.  The 
geomaterials  of  concern  here  are  generally  partially  satmated.  Volumetric  strain  for  these 
materials  is  determined  by  combining  measurements  of  the  sample  deformations  based  on 
an  assumed  shape.  The  method  used  to  calculate  volumetric  strain  must  be  understood  so 
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that  consistency  is  maintained  between  the  mathematical  model  and  the  method  used  to 
characterize  the  material  property.  Here  and  throu^out  this  chapter,  Lagrangian,  or  total, 
notation  is  used  to  define  strain,  i.e.,  strain  is  equal  to  the  change  in  dimension  divided  by 
the  original  dimension,  and  compression  is  taken  to  be  positive. 

Most  mathematical  models  assume  small  strain,  at  least  in  the  increment,  and 
therefore  the  most  common  method  for  calculating  volumetric  strain  is 

8=e+e+s  32 

V  X  y  Z 

where  is  the  volumetric  strain  and  e^,  e^,  and  are  the  axial  strains  in  the  three 
orthogonal  directions  assuming  cartesian  coordinates.  In  cylindrical  coordinates  the 
individual  strains  are  89,  and  e„  with  8^  and  89  often  assumed  to  be  equal. 

Other  methods  for  calculating  volumetric  strain  based  on  an  assumed  shape  were 
presented  by  Ehrgott  (1971)  for  cylindrical  specimens.  The  shape  of  the  specimen  during 
hydrostatic  compression  is  influenced  by  the  end  restraint.  If  the  specimen  does  not  move 
fi'eely  at  the  ends,  then  the  shape  is  similar  to  a  double  truncated  cone  and  the  volumetric 
strain  is  calculated  as 


If  the  specimen  ends  are  free  to  move,  the  shape  is  more  like  that  of  a  cylinder  and  the 
volumetric  strain  is  calculated  as 

8^  -  8^  +  2  8^  -  2  8^  8^  +  8^  (8^  -  1)  3.4 

If  the  higher  order  terms  of  Equation  3.4  are  removed,  it  reduces  to  Equation  3.2  for  small 
strains. 
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Given  a  right  rectangular  shape  as  in  Figure  3.6,  the  original  volume  can  be 
written  as 


^  L  RP 

O 

and  the  current  volume  can  be  defined  as 

V.  =  (L  -  AL)  (R  -  AR)  (P  -  AP) 

=  LRP(l  (I  -  Bp) 


3.5 


3.6 


where  L,  R  and  P  are  the  lengths  of  each  side,  A  is  the  change  in  dimension  and  e  is  the 
total  strain.  Defining  as 


e 


V 


=  1 


3.7 


and  substituting  Equations  3.5  and  3.6  gives 


®v  =  1  - 


Z  i?  P  (1  -  (1  -  e^)  (1  -  ep 

L  RP 


3.8 


Equation  3.8  is  similar  to  Equation  3.4  and,  if  the  higher  order  terms  are  removed,  is  the 
same  as  the  small  strain  definition  in  Equation  3.2.  Figure  3.7  shows  a  calculation  of 
volumetric  strain  using  Equation  3.2  for  the  small  strain  assumption  and  Equation  3.8  for 
the  large  strain  assumption,  and  assuming  isotropic  deformation  such  that  the  axial  strains 
are  equal.  Most  brittle  geomaterials  will  exhibit  a  volumetric  strain  during  hydrostatic 
compression  of  less  than  10  percent,  and  certainly  less  than  20  percent.  Figure  3.7  shows 
that  the  using  the  small  strain  assumption  is  only  3  percent  larger  than  the  8^  using  the 
large  strain  assumption  at  an  axial  strain  of  about  0.03,  and  only  7  percent  larger  at  an 
axial  strain  of  about  0.07.  These  axial  strains  correspond  to  volumetric  strains  of  about 
0.1  and  0.2  (10  and  20  percent),  respectively. 


Equation  3.2  (small  strain  assumption) 
Equation  3.8  (large  stran  assumption) 
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Brittle  geomaterials  exhibit  nonlinear  stress-strain  behavior  during  application  of 
hydrostatic  pressure  as  shown  in  Figure  3.8,  where  volumetric  strain  is  calculated  using 
Equation  3.2  and  the  pressure  is  the  true  stress  being  applied  to  the  concrete  specimen.  In 
Figure  3.8.a  shows  the  results  from  hydrostatic  compression  tests  on  four  concrete 
specimens.  The  initial  compression  slope  of  the  hydrostatic  response  is  fairly  constant 
and  the  material  is  essentially  elastic.  The  response  then  becomes  concave  to  the  strain 
axis,  but  then  reverses  to  become  concave  to  the  stress  axis.  This  behavior  indicates  a 
change  is  occurring  within  the  material,  but  it  is  not  failing,  i.e.,  when  pressure  is 
removed  the  material  retains  a  shape  with  no  easily  visible  cracks  or  damage.  During 
unloading,  the  slope  is  nearly  constant  initially,  but  becomes  concave  to  the  stress  axis  as 
pressure  approaches  zero,  which  is  another  indication  that  the  material  has  experienced 
some  internal  change.  Figure  3.8.b  shows  the  results  from  hydrostatic  tension  tests  on 
two  concrete  samples  and  hydrostatic  compression  tests  on  three  samples.  In  tension,  the 
initial  slope  is  a  linear  extension  of  the  initial  compressive  slope.  As  the  tensile  stress 
increases,  the  curve  becomes  more  concave  to  the  strain  axis  and  the  material  begins  to 
experience  more  internal  changes.  In  tension,  an  easily  visible  failure  occurs  when  the 
material  separates.  Little  information  is  available  to  describe  “unloading”  during 
hydrostatic  tension  because  the  failure  of  the  material  is  sudden  and  catastrophic  and 
control  of  the  test  is  often  lost. 

The  idealized  brittle  geomaterials  here  are  assumed  to  consist  of  relatively  hard 
granular  particles  cemented  together  by  a  weaker  paste,  whether  introduced  by  man  as 
with  concrete  or  by  nature  as  in  cemented  sands  and  rocks,  with  other  inclusions  such  as 
micro-  and  macrocracks  and  voids.  The  voids  may  be  filled  with  liquid  such  as  water  or 
gas  such  as  air.  During  application  of  the  hydrostatic  pressure  the  brittle  geomaterial 
specimen  undergoes  several  physical  changes.  These  changes  are  illustrated  in  Figure  3.9 
where  compression  is  positive.  In  region  A  of  Figure  3.9,  the  material  is  essentially 
elastic  and  the  stresses  are  not  sufficient  to  cause  slippage,  significant  coalescence  of 
cracks,  or  significant  damage  to  the  paste  or  particles.  The  pressure  being  applied  to  the 
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Figure  3.8.  Response  of  concrete  to  hydrostatic  pressure. 
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Figure  3.9.  Material  changes  during  application  of  hydrostatic  stress. 
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specimen  is  easily  carried  by  the  material  constituents.  In  region  B,  crack  formation  and 
coalescence  of  cracks  within  the  paste  begins  to  occur.  The  cracks  formed  and  extended 
in  region  B,  along  with  the  existing  voids,  will  allow  for  the  sliding  and  particle 
rearrangement  that  occurs  in  region  C  and  additional  breakdown  of  the  paste.  As  the 
pressure  increases  in  region  C,  some  void  and  crack  closure  is  possible  and  cracking  and 
breakdown  of  the  granular  particles  may  begin.  In  region  D,  as  pressure  increases,  less 
particle  rearrangement  occurs  while  void  and  crack  closure  becomes  significant.  At  the 
extreme  high  pressures  in  region  D,  the  voids  and  openings  from  the  cracks  are  closed  or 
filled  with  fluid  such  that  the  applied  pressures  are  being  carried  by  the  pore  fluid  and  not 
the  skeleton  of  the  material. 

During  the  application  of  all  around  tensile  stresses,  the  response  of  the  material, 
regions  E  and  F  in  Figure  3.9,  can  be  dramatic  in  that,  unlike  compression,  catastrophic 
failure  of  the  test  specimen  may  occur.  Region  E  is  similar  to  region  B  in  compression 
except  that  the  majority  of  the  changes  to  the  material  are  due  to  crack  propagation  and 
coalescence.  At  some  level  of  tension,  significant  cracks  will  form  and  the  material  will 
fail  with  rapid,  but  not  immediate,  loss  of  capacity  (region  F).  The  loss  is  not  immediate 
because  the  loss  of  strength  will  begin  prior  to  complete  separation  at  the  significant 
crack.  Since  complete  separation  ultimately  occurs,  the  specimen  cannot  retain  a  residual 
strength. 

3.3  DEVIATORIC  LOADING 

A  deviatoric,  or  shear,  loading  is  often  applied  to  a  test  specimen  after  some  level 
of  compressive  (HC)  or  tensile  hydrostatic  (HT)  pressure  has  been  applied.  The  load,  or 
deformation,  is  changed  while  the  response  of  the  material  is  monitored.  Desai  and 
Siriwardane  (1984)  describe  several  of  the  different  types  of  loading  paths  that  might  be 
applied.  The  states  of  principal  stress  for  triaxial  experiments  on  cubic  and  cylindrical 
specimens  are  shown  in  Figure  3.10.  In  experiments  on  cylindrical  specimens  the  stresses 
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Oj  and  O3  are  assumed  to  be  equal.  Stress-paths  for  typical  mechanical  property 
experiments  are  illustrated  in  the  triaxial  plane  in  Figure  3.11.  In  most  triaxial  tests  the 
specimen  is  subjected  to  an  initial  confining  stress  such  that  the  principal  stresses  are  all 
equal,  0^  =  0^  =  02  =  Oy  In  the  conventional  triaxial  compression  (CTC)  test  two  of  the 
stresses,  usually  Oj  and  O3,  are  held  constant  while  Oj  is  increased  until  failure  of  the 
specimen  occurs.  The  deviator  stress  and  hydrostatic  pressure  both  increase  during  the 
loading.  In  the  reduced  triaxial  compression  (RTC)  test,  Oj  and  O3  are  reduced  while  o,  is 
held  constant  so  that  the  shear  stress  increases  while  the  pressure  decreases.  In  the 
triaxial  compression  (TC)  test,  the  stress  is  applied  so  that  the  pressure  remains  constant 
while  the  shear  stress  increases.  Similar  loadings  can  be  applied  so  that  the  shear  stress 
decreases  to  conduct  experiments  in  extension  (RTE,  CTE,  and  TE).  In  the  simple  shear 
(SS)  test,  the  stresses  are  applied  such  that  the  mean  pressure  is  kept  constant.  In  a 
proportional  loading  (PL)  test,  the  loading  of  the  specimen  is  conducted  while 
maintaining  a  constant  ratio  between  o„  Oj  and  to  give  a  particular  loading  path. 

The  stress  reported  from  experiments  involving  deviatoric  loads  can  either  be  the 
true  stress  or  the  total  stress  (conjugate  stress  to  strain  based  on  the  original  dimensions). 
The  total  stress  may  be  provided  when  some  deformation  measurements  are  not  recorded. 
On  cylindrical  specimens,  if  the  lateral  deformations  are  not  recorded  during  an 
unconfined  compression  test,  the  load  applied  to  the  specimen  may  be  divided  by  the 
original  cross-sectional  area  to  give  the  total  stress  for  the  material  response.  In  tests 
where  pressure  (a  true  stress)  is  applied  to  the  specimen,  the  lateral  deformations  are 
generally  recorded  and  the  axial  stress  is  calculated  based  on  the  current  dimensions  of 
the  specimen  to  give  a  true  stress.  Throughout  this  chapter,  stresses  are  assumed  to  be  the 
true  stresses  based  on  the  current  dimensions  of  the  specimen  and  strains  are  the 
engineering  strains  based  on  the  original  dimensions.  The  use  of  engineering  strain 
allows  the  material  response  to  be  presented  as  a  function  of  the  total  strain  of  the 
specimen. 
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The  amount  of  data  on  the  response  of  brittle  geomaterials  to  the  pressures 
experienced  during  projectile  penetration  (see  Figures  3.3  through  3.5)  is  very  limited  and 
most  of  the  available  data  is  in  compression.  The  majority  of  these  tests  are  conducted  as 
conventional  triaxial  compression  because  they  provide  a  straightforward  determination 
of  the  elastic  constants  (shear  modulus  G,  Young’s  modulus  E,  and  Poisson’s  ratio  v)  and 
they  are  easy  to  perform. 

3.3.1  Compression 

Chen  (1982)  gives  a  description  for  the  typical  response  of  a  concrete  to  uniaxial 
compression.  A  typical  stress-strain  curve  for  CSPC  concrete  in  uniaxial  compression  is 
shown  in  Figure  3.12.  The  response  is  nearly  linear  elastic  initially  since  the  cracks  and 
voids  existing  in  the  specimen  are  nearly  unchanged.  Between  a  stress  of  about  30 
percent  and  75  to  90  percent  of  the  maximum  compressive  strength  the  response  curves 
gradually  increase  in  curvature.  Above  this  region  the  response  curves  bend  sharply  and 
the  specimen  approaches  its  peak  strength.  Bond  cracks  caused  by  tensile  stress 
concentrations  between  binder  and  aggregate  (van  Vliet  and  van  Mier  1996)  start  to 
extend  at  stresses  between  30  to  50  percent  of  the  peak.  The  crack  propagation  is  stable 
in  that  the  crack  lengths  quickly  reach  their  final  values  if  the  applied  stress  is  kept 
constant.  At  stresses  between  50  to  75  percent  of  the  peak  some  cracks  begin  to  bridge 
into  the  binder  while  some  bond  cracks  continue  to  grow.  Above  about  75  to  90  percent 
of  the  peak  stress,  the  largest  cracks  reach  their  critical  lengths  and  the  available  internal 
energy  of  the  system  is  greater  than  the  required  crack-release  energy,  or  the  energy 
required  to  propagate  a  crack.  The  system  is  unstable  and  complete  failure  can  occur 
even  if  the  load  is  kept  constant.  In  Figure  3.12  at  this  point  the  lateral  deformation  is 
expanding  more  rapidly  than  the  axial  compression  and  the  specimen  begins  to  dilate.  A 
posttest  description  of  the  test  specimen  describes  a  diagonal  shear  across  the  specimen 
firom  top  to  bottom. 
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The  softening,  or  descending,  portion  of  the  stress-strain  curve  is  greatly 
influenced  by  specimen  size  and  boundary  conditions.  Softening  occurs  when  the 
microcracks  that  began  forming  during  the  pre-peak  portion  of  the  experiment  coalesce  to 
form  a  damage  zone  that  greatly  weakens  the  specimen  (Jansen  and  Shah,  1995).  Since 
localization  occurs  within  the  damage  zone,  the  post-peak  portion  of  the  stress-strain 
curve  is  a  material  property,  but  dependent  on  the  specimen  length  over  which  the 
deformation  is  measured.  Van  Vleit  and  van  Mier  (1996)  conclude  that  softening  is  a 
combined  material  and  structural  property.  If  the  length  is  confined  to  the  area  of  the 
localization,  the  effect  of  specimen  length-to-diameter  ratio  on  the  post-peak  energy  is 
removed  for  ratios  between  2.0  to  5.5.  Results  from  the  experiments  by  Janzen  and  Shah 
showed  that  the  fracture  zone  is  200  mm  or  more  long  for  100  mm  diameter  specimens. 
Figure  3.13  shows  the  stress-strain  behavior  for  different  length-to-diameter  ratio 
specimens  for  a  normal  concrete  with  an  average  UC  strength  of  about  48  MPa.  The 
ascending  portion  of  the  curves  are  the  same  for  all  ratios.  The  descending  portion  of  the 
curve  is  steeper  for  the  high  ratio  specimens  since  the  strain  is  based  on  the  original 
length  of  the  specimen,  but  the  deformation  is  occurring  within  the  200  mm  long  damage 
zone.  During  softening,  the  cracks  and  shear  bands  will  grow  to  a  size  comparable  to  the 
characteristic  dimension  of  the  specimen  (van  Vliet  and  van  Mier,  1996). 

Experiments  by  van  Vliet  and  van  Mier  (1996)  also  showed  that  a  specimen 
length  at  least  twice  the  lateral  dimension  is  required  to  avoid  size  effects  in  peak  strength 
with  high  fiiction  end  restraints,  and  that  the  post-peak  deformation  is  independent  of  size 
for  ratios  from  0.5  to  2.0.  A  length  to  width  ratio  of  at  least  2.0  should  be  used  to  avoid 
influencing  the  pre-peak  and  ultimate  strength  of  the  material.  Strain  calculated  in  the 
softening  branch  of  the  response  should  be  based  on  the  characteristic  length  for  the 
specimen,  which  ideally  is  the  same  as  the  specimen  height.  In  order  to  capture  the  pre¬ 
peak  and  post-peak  response  and  the  peak  strength,  a  length  to  lateral  dimension  of  two 
should  be  used  for  conventional  concrete. 
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Additional  experimental  results  by  Lee  and  Wiliam  (1995)  further  illustrate  the 
effects  of  specimen  size.  They  conducted  uniaxial  compression  experiments  on 
cylindrical  concrete  specimens  with  aspect  ratios  of  1 .8, 1 .2,  and  0.6  under  deformation 
control.  The  ultimate  strength  of  the  1.8  and  1.2  aspect  ratio  specimens  was  the  same  and 
the  strength  of  the  0.6  aspect  ratio  specimen  was  about  30  percent  greater.  They 
measured  both  axial  and  radial  deformation  and  note  that  the  initial  Poisson’s  ratio  of 
0.18  in  the  prq)eak  region  increased  to  about  12  (based  on  deformation)  in  the  post-peak 
regions.  This  “fictitious  Poisson’s  effect”  was  the  same  for  all  specimen  heights  and  was 
constant  in  the  post-peak  regime.  The  post-peak  response  is  the  result  of  localized 
splitting  that  initiates  at  the  peak  strength,  and  the  energy  dissipation  is  the  combined 
effect  of  axial  splitting  and  shear  dissipation.  The  characteristic  length  for  the  test 
specimens  was  estimated  to  be  0.095  h  based  on  fracture  energies  for  mode  I  and  mode  II 
type  failures  and  0.27  h  to  0.44  h  based  on  analytic  solutions  of  fracture  mechanics,  where 
h  is  the  specimen  height. 

Similar  observations  have  been  made  for  rock  (Bruno  and  Nelson,  1991;  Desai,  et 
al,  1990;  Sulem  and  Vardoulakis,  C&D).  Large  rock  masses  contain  faults,  joints,  and 
changes  in  lithology  across  bedding  planes  that  separate  the  different  strata.  Individual 
specimens  of  rock  are  a  complex  assembly  of  crystals,  grains,  smaller  matrix  materials, 
cementing  material,  and  porespace  (Bruno  and  Nelson,  1991).  The  observed  inelastic 
deformations  in  material  property  experiments  are  the  result  of  microstructural  damage 
and  fracturing  and  are  similar  to  the  results  from  experiments  on  concrete. 

3.3.2  Tension 

The  tensile  response  of  a  material  generally  provides  the  basic  input  for  fracture 
mechanics  models  (Hordijk,  et.al.  1989).  Like  in  compression,  the  crack  evolution 
caimot  be  obtained  directly,  but  must  be  derived  from  the  complete  stress-displacement 
relation.  The  ultimate  tensile  strength  is  determined  by  the  strength  and  stiffiiess  of  the 
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particles  and  the  cementing  material  and  the  bond  strength  between  the  two.  In  cases 
when  the  particles  and  cement  are  similar  in  strength,  the  crack  paths  pass  through  the 
aggregate  particles.  This  results  in  a  more  brittle  response  than  if  the  crack  path  is  more 
tortuous  and  must  travel  around  particles  instead  of  through  them. 

Rossi,  et  al  (1994)  present  results  from  direct  tensile  experiments  on  cylindrical 
specimens  of  concrete  with  different  diameters  but  all  with  a  height-to-diameter  ratio  of 
two.  They  note  a  size  effect  in  that  the  strength  for  the  smaller  specimens  is  higher  than 
that  for  the  larger  specimens.  Specimens  of  the  same  size  also  exhibited  a  Yoxmg’s 
modulus  that  was  equal  to  the  modulus  in  compression.  The  scale  effect  is  hypothesized 
to  be  a  function  of  the  specimen  volume  and  the  total  volume  of  the  coarsest  aggregate. 
The  basis  is  that  the  weakest  link  in  the  material  is  the  cement  paste  that  contains  faults 
such  as  bubbles,  microcracks,  initial  stresses,  etc.,  so  that  crack  initiation  begins  in  the 
paste. 


Results  from  a  direct  pull  tension  experiment  on  a  concrete  specimen  with 
height/diameter  ratio  of  about  2  are  shown  in  Figure  3.14.  The  experiment  was 
conducted  imder  load  control  therefore  the  catastrophic  failure  of  the  specimen  resulted  in 
no  post-peak  data.  In  uniaxial  tension,  bond  microcracks  begin  to  grow  in  concrete  at 
stresses  above  about  60  percent  of  the  ultimate  tensile  strength.  The  onset  of  unstable 
crack  propagation  occurs  at  about  75  percent  of  the  ultimate  strength  and  the  stress-strain 
response  becomes  more  concave.  The  crack  direction  is  generally  perpendicular  to  the 
direction  of  loading.  The  ultimate  strength  of  concrete  m  uniaxial  tension  is  about  5  to  10 
percent  of  the  strength  in  compression  (Chen,  1982). 

3.3.3  Effect  of  Confinement 

Confinement  has  a  significant  effect  on  the  response  of  brittle  geomaterials  to 
deviatoric  loads.  Figure  3.15  schematically  shows  the  gross  effect  of  confinement  on 
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material  response  in  compression.  As  discussed  above,  at  low  confining  pressures  these 
materials  behave  in  a  brittle  manner  where  the  propagation  of  unstable  cracks  leads  to 
their  coalescence  and  the  eventual  formation  of  failure  planes  (Ladanyi  and  Aubertin, 
1990;  Lubada,  et.al.,  1996;  Ashby  and  Hallam,  1986).  The  failure  plane  formation  is 
illustrated  by  dilation  in  the  overall  response  of  the  material.  As  confining  pressure 
increases,  the  material  response  to  shear  loading  transitions  through  a  brittle  to  ductile 
regime  where  dislocation  motion  and  fracture  propagation  processes  interact  and  the 
formation  of  a  shear  band  and  strain  softening  are  suppressed.  Microcracks  still  occm 
and  produce  dilation,  but  ductile  processes  at  the  crystalline  scale  are  also  influencing  the 
material  response  (Ladanyi  and  Aubertin,  1990).  Increasing  confining  pressure  results  in 
increasing  strength.  At  high  confining  pressures  deformation  occurs  at  the  intercrystalline 
level  resulting  in  a  ductile  response.  The  effect  of  confining  pressure  on  the  TC  response 
of  concrete  specimens  is  shown  in  Figure  3.16.  As  confining  pressure  increases,  the 
response  transitions  fi-om  a  brittle  response  at  low  confining  pressure  to  a  ductile  response 
at  high  confining  pressure.  Posttest  descriptions  of  the  test  specimens  are  similar  to  those 
in  Figure  3.15.b  and  c. 

As  discussed  in  Section  3.2  and  illustrated  in  Figure  3.9,  as  pressure  is  applied  to 
brittle  cemented  geomaterials  it  is  hypothesized  that  the  material  is  altered  through  the 
creation  and  propagation  of  dislocations,  rearrangement  of  particle  orientation,  and 
compaction.  To  test  this  hypothesis,  TC  experiments  were  conducted  on  concrete 
specimens  at  two  different  confining  pressures  where  some  of  the  specimens  were 
hydrostatically  pre-stressed  prior  to  the  TC  experiment.  Results  from  TC  experiments  at 
confining  pressures  of  25  and  100  MPa  are  shown  in  Figures  3.17  and  3.18,  respectively. 
For  the  25  MPa  experiments,  one  pre-stress  was  conducted  to  100  MPa,  or  at  about  the 
knee  of  the  HC  pressure- volume  response  where  breakdown  and  cracking  in  the  matrix  is 
beginning,  and  one  pre-stress  was  conducted  to  400  MPa,  or  at  a  level  of  pressure  where 
the  material  has  undergone  significant  changes.  The  shear  response  of  the  specimen  pre¬ 
stressed  to  100  MPa  is  similar  to  the  response  of  the  specimen  with  no  pre-stress,  but 


Figure  3.16.  Effect  of  confining  pressure  on  results  from  TC  tests  on  concrete  specimens, 
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does  show  a  tendency  toward  dilation  earlier  in  the  shear  loading.  The  shear  response  of 
the  specimen  pre-stressed  to  400  MPa  is  much  softer  than  the  specimen  with  no  pre¬ 
stress,  has  a  strength  slightly  lower,  and  shows  significant  dilation  during  shear  that 
begins  almost  as  soon  as  the  shear  loading  is  applied.  The  appearance  of  the  specimens 
after  the  TC  experiment  was  similar  for  all  pre-stresses.  For  the  100  MPa  experiments 
(Figure  3.18),  a  pre-stress  to  400  MPa  was  applied  to  two  specimens.  The  shear  response 
of  the  pre-stressed  specimens  is  stiffer  than  the  specimens  with  no  pre-stress  and  has  an 
ultimate  strength  that  is  slightly  higher.  The  pre-stressed  specimens  begin  to  dilate  soon 
after  the  shear  loading  is  applied.  Test  specimens  that  had  been  pre-stressed  to  400  MPa 
exhibited  slightly  more  bulging  near  the  midheight  and  a  more  visible  shear  plane  after 
the  TC  experiment. 

3.3.4  Effect  of  Strain  Rate 

The  dynamic  response  of  the  target  materials  is  of  interest  since  the  stress  field 
within  the  target  during  the  penetration  process  is  being  applied  over  a  period  of  only  a 
few  milliseconds.  Conventional  testing  equipment  using  pressure  chambers  and  loading 
rams  can  be  used  to  test  materials  at  strain  rates  to  about  10  /s.  These  devices  typically 
use  fast  opening  valves  or  explosives  to  drive  rams  that  apply  the  stresses  to  the  test 
specimen.  Measurements  are  made  using  mechanical  based  systems  such  as  strain  gages, 
linear  variable  differential  transformers  (LVDTs),  coiled  wire  pressure  gages,  etc.  Drop 
hammers  and  anvils  can  also  be  used  to  conduct  experiments  to  strain  rates  of  about  10  /s. 
Above  this  rate  of  loading  wave  propagation  effects  between  the  specimen  and  test  device 
must  be  taken  into  account  (Sierakowski  1984).  Pressure  bars  such  as  the  split 
Hopkinson  bar  have  been  used  to  obtain  data  at  strain  rates  up  to  about  10^  /s.  Data  from 
the  experiments  must  be  analyzed  using  wave  propagation  theories  in  order  to  deduce  the 
stress  and  strain  within  the  test  specimen.  Another  wave  propagation  experiment  that  can 
provide  data  at  strain  rates  greater  than  10^  /s  is  the  flyer  plate  test. 
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The  test  methods  have  been  applied  by  many  researchers  to  investigate  the 
dynamic  response  of  brittle  geomaterials  to  compressive  and  tensile  states  of  stress. 
Results  from  dynamic  compression  experiments  conducted  by  several  researchers  were 
summarized  by  Bischoff  and  Perry  (1986).  The  ratio  of  d5mamic  to  static  compressive 
strength  is  shown  versus  the  log  of  strain  rate  in  Figure  3.19.  Quasi-static  loading  is 
generally  taken  to  be  about  10'^  to  lO"^  /s.  The  data  show  a  gradual  increase  in  strength  to 
a  strain  rate  between  10"'  and  1  /s,  where  the  dynamic  strength  begins  to  increase  rapidly 
as  the  strain  rate  increases  to  a  ratio  of  almost  2  at  a  rate  of  10^  /s.  Hughes,  et.al.  (1993) 
show  the  effects  of  strain  rate  on  the  tensile  strength  of  concrete  in  Figure  3.20. 
Comparison  of  Figures  3.19  and  3.20  indicates  that  the  increase  in  strength  with  strain 
rate  is  more  dramatic  in  tension  than  in  compression. 

Confinement  of  the  test  specimen  has  been  shown  to  influence  the  increase  in 
strength  as  strain  rate  increases  (Gran,  et.  al.  1989).  Results  by  Gran,  et.  al.  (1989)  in 
Figure  3. 21. a  for  a  concrete  with  static  imconfined  compressive  strength  of  about  100 
MPa  indicate  that  the  ratio  between  the  static  and  dynamic  strength  (1.3  to  1.4)  is  nearly 
constant  as  pressure  increases  for  dynamic  strain  rates  of  about  1.3  /s  to  5  /s.  Results  by 
Yamaguchi,  et.  al.  (1989)  in  Figure  3.21.b  for  a  concrete  with  unconfined  compressive 
strength  of  about  20  MPa  indicate  that  the  dynamic  strength  increases  gradually  from  the 
static  strength  as  pressure  increases  for  strain  rates  of  0.254  /s  and  lower.  Results  from 
static  and  dynamic  triaxial  compression  experiments  conducted  at  WES  on  CSPC 
concrete  are  shown  in  Figure  3.22  where  strain  rate  (approximately  0.1  /s)  is  the  axial 
strain  at  peak  stress  divided  by  the  time  to  peak,  all  adjusted  to  the  start  of  the  shear  phase 
of  the  experiment.  The  data  show  the  ratio  between  dynamic  and  static  strength 
decreasing  slightly  as  the  confining  pressure  increases,  but,  at  this  dynamic  strain  rate,  the 
data  are  within  the  scatter  of  the  data  in  Figure  3.19. 

Loading  rate  also  affects  the  shape  of  the  stress-strain  response  of  brittle 
geomaterials.  Figure  3.23  shows  the  effect  of  strain  rate  on  the  unconfined  compression 
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Figure  3.20.  Strain  rate  effects  on  concrete  tensile  strength  (Hughes,  et.  al  1993). 
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Figure  3.21 .  Effect  of  pressure  on  dynamic  strength  of  concrete. 
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stress-strain  response  of  three  concretes.  Figure  3.24  shows  the  effect  of  strain  rate  on  the 
stress-strain  response  of  CSPC  concrete  at  confining  pressures  of  about  0, 10,  and  20 
MPa.  The  pre-peak  portion  of  the  dynamic  responses  are  generally  stiffer  than  the  static 
response,  but  the  post-peak  softening  response  tends  to  become  softer  (more  gradual)  as 
the  strain  rate  increases.  The  strain  at  peak  strength  does  not  change  significantly  as  the 
strain  rate  increases.  Similar  observations  can  be  made  for  the  tensile  response  of 
concrete  (Figure  3.25). 

It  has  been  suggested  that  the  enhanced  performance  of  brittle  geomaterials  during 
dynamic  loading  is  due  to  lateral  inertial  confinement  of  the  interior  of  the  specimen. 
Tang,  et.  al.  (1992)  conducted  analyses  that  showed  the  magnitude  of  this  confinement  is 
only  about  0.1  to  0.2  percent  of  the  applied  axial  stress  and  could  not  account  for  the 
significant  strength  increase.  The  development  of  internal  microcracks  is  responsible  for 
the  increasing  non-linearity  of  brittle  geomaterials  when  loaded  quasi-statically.  The 
resistance  of  these  cracks  to  propagate  during  increased  loading  rate  is  believed  to  be  die 
cause  of  the  increased  strength  (Bischoff  and  Perry  1986).  During  high  loading  rates  the 
internal  cracks  do  not  have  time  to  follow  the  path  of  least  resistance  and  are  forced  to 
propagate  along  paths  that  go  through  aggregate  and  paste  that  would  normally  be 
avoided.  Higher  strength  concretes  do  not  show  as  significant  an  increase  in  strength  as 
loading  rate  increases  since  cracks  are  often  forced  to  propagate  through  aggregate  even 
during  quasi-static  loading.  As  shown  in  Figures  3.19  and  3.20,  a  wide  scatter  exists  in 
the  data.  The  discrepancies  are  caused  by  several  factors  including  the  inherent 
heterogeneity  of  concrete,  the  concrete  mix  proportions,  aggregate  type,  curing 
conditions,  age  at  testing,  specimen  geometry,  and  boimdary  conditions  during  the 
experiment  (Fu,  et.  al.  1990).  Many  of  these  same  observations  can  also  be  applied  to 
other  brittle  geomaterials  such  as  rock. 
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Figure  3.24.  Effect  of  confining  pressure  on  the  static  and  dynamic  stress-strain  response 
of  CSPC  concrete. 
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Figure  3.25.  Rate  effect  on  load-deformation  curve  in  tension  (Kormeling,  et.  al.  1987  in  Weerheijm,  et.  al.  1993). 
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CHAPTER 4 

NONLINEAR,  INELASTIC  FRACTURE  MODEL 

4.1  MODELING  OF  BRITTLE  GEOMATERIALS 

Two  procedures  are  generally  employed  to  formulate  constitutive  theories  for  the 
nonlinear  stress-strain  response  of  materials  (Rohani  and  Thompson  1970).  One 
procedure  seeks  to  find  a  physical  phenomaia  that  is  the  outcome  of  a  framework  of 
elaborate  mathematics.  Although  the  mathematics  is  complete,  it  may  be  difficult  to  find 
actual  materials  whose  response  is  simulated  by  the  predicted  stress-strain  curves.  The 
other  procedure,  often  referred  to  as  the  physical  approach,  seeks  to  find  the  appropriate 
mathematical  forms  needed  to  model  the  observed  physical  behavior  of  actual  materials. 
The  essential  material  responses  are  incorporated  into  the  theory  at  the  beginning.  One 
pleasing  feature  of  this  approach  is  that  new  criteria  or  mathematical  forms  can  be 
incorporated  into  the  theory  as  appropriate.  Since  no  one  approach  or  theory  can  be 
expected  to  describe  the  response  of  a  material  under  all  conditions,  the  approach  taken 
may  come  down  to  the  preferences  of  the  model  developer.  Also,  the  desirable  features 
of  both  approaches  might  be  used  in  the  model  development. 

Desai  and  Siriwardane  (1984)  provide  seven  axioms  that  the  constitutive  model 
must  obey.  The  axiom  of  determinism  states  that  current  state  of  stress  and  strain  due  to 
an  external  force  is  dependent  on  the  history  of  forces  that  have  been  experienced.  The 
axiom  of  causality  states  that  a  material  response  cannot  occur  without  a  cause.  The 
axiom  of  objectivity  states  that  the  material  properties  cannot  vary  with  the  motion  of  the 
observer.  The  axiom  of  neighborhood  states  that  the  response  at  a  point  is  not  affected  by 
conditions  that  are  far  away.  The  axiom  of  memory  states  that  the  constitutive  variables 
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are  not  affected  by  the  values  of  the  constitutive  variables  in  the  distant  past.  The  axiom 
of  equipresence  states  that  the  variables  identified  for  a  material  should  be  present  in  all 
of  the  constitutive  equations  for  that  material.  The  axiom  of  admissibility  states  that  the 
constitutive  equations  must  be  consistent  with  the  physical  laws,  such  as  continuity, 
motion,  momentum,  entropy,  etc.  Since  a  constitutive  model  is  intended  to  relate  a 
physical  response  to  applied  loads,  the  model  must  obey  the  above  axioms  that  govern  the 
response. 

The  various  approaches  used  to  model  geomaterials  can  be  foxmd  in  several  texts 
(Chen  1982,  Chen  and  Saleeb  1982,  Desai  and  Siriwardane  1984,  Chen  and  Baladi  1985). 
The  approaches  include  linear  and  nonlinear  elasticity  with  and  without  failure  criteria 
and  with  and  without  fi'acture,  hypoelasticity,  hyperelasticity,  and  plasticity  with  and 
without  fracture,  as  well  as  endochronic,  viscoelastic,  and  viscoplastic  methods.  The 
linear  elastic  models  are  widely  used  in  practice  because  tiiey  are  relatively  easy  to  use 
and  understand  and  the  model  constants  are  usually  easy  to  obtain.  The  stress-strain 
relationships  for  most  geomaterials  however  can  be  better  modeled  by  using  nonlinear 
elasticity  equations.  The  elastic  methods  work  well  when  the  loading  is  monotonic  and 
well  below  failure.  If  the  effects  of  unloading  and  reloading  are  significant  then  a 
separate  criteria  for  these  conditions  may  be  included.  The  plasticity  based  methods  use 
yield  surfaces  and  flow  criteria  to  determine  the  material  response  incrementally.  They 
cm  range  fi’om  the  simple  elastic-perfectly  plastic  model  to  complicated  with  multiple 
response  surfaces  and  flow  criteria. 

In  recent  years,  efforts  to  use  fi-acture  and  damage  mechanics  methods,  both  alone 
and  in  conjimction  with  elasticity  and  plasticity,  to  describe  the  response  of  brittle 
geomaterials  have  increased  (Elfgren  1989).  Fracture  mechamcs  evolved  as  a  result  of 
failures  in  metal  structures  built  in  the  19*  century  through  World  War  II  (Ewalds  and 
Wanhill  1984).  Investigations  revealed  that  material  deficiencies  could  initiate  cracking 
and  fi'acture  at  flaws  and  stress  concentrations.  Similarly,  brittle  geomaterials  contain 
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microcracks,  pores,  voids,  and  stress  concentrations  where  fractures  can  initiate.  The  pre¬ 
failure  stress-strain  behavior  is  often  characterized  by  models  based  on  the  elasticity  and 
plasticity  theories  of  continuum  mechanics  (Y amaguchi  and  Chen  1987),  but  the  tensile 
and  compressive  failures  are  primarily  due  to  macroscopic  crack  propagation  (Jenq  and 
Shah  1987)  which  cannot  be  explained  by  using  only  these  theories.  Fracture  mechanics 
is  concerned  almost  entirely  with  the  fracture-dominant  behavior  characterized  by  highly 
localized  plasticity  and  essentially  macroscopic  sized  defects  (Ewalds  and  Wanhill  1984). 

A  damage  model  is  one  that  contains  a  specific  internal  variable  that  qualifies  the 
state  of  damage  locally  and  a  kinetic  equation  that  defines  the  evolution  of  the  damage 
with  the  applied  load  or  time  (Krajcinovic  1984).  The  internal  variable  is  tracked  during 
loading  and  used  to  calculate  or  trigger  the  calculation  of  a  damage  variable.  The  damage 
variable  usually  varies  from  zero  to  one  and  is  applied  to  the  stifBiess  of  the  material. 

The  damage  variable  can  be  either  a  scalar  for  inducing  isotropic  damage  or  a  vector  for 
inducing  anisotropic  damage.  As  the  damage  is  accumulated,  the  stiffiiess  of  the  material 
is  softened  to  effectively  induce  a  strain-softening  effect.  Individual  models  can  vary  on 
the  internal  variable  that  is  tracked  and  how  the  damage  is  calculated  and  applied. 

4.2  NIF  MODEL  DEVELOPMENT 

A  non-linear,  inelastic  fracture  (NIF)  model  will  be  developed  to  capture  some  of 
the  response  characteristics  of  brittle  geomaterials  to  the  severe  loading  conditions  that 
occur  during  high-velocity  projectile  impact.  In  this  sense,  the  model  will  be  developed 
following  the  physical  approach  to  replicate  the  responses  presented  in  the  previous 
chapter  assuming  a  frue  stress-total  strain  representation.  The  state  of  stress  in  the  model 
will  be  determined  uniquely  by  the  current  state  of  strain.  The  model  will  be  nonlinear 
and  inelastic,  include  the  tensile  characteristics  of  the  material,  and  capture  the  brittle 
shear  response  at  low  pressures,  the  ductile  shear  response  at  high  pressures  and  the 
transition  between  the  brittle  and  ductile  behavior. 
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The  behavior  of  a  material  is  modeled  mathematically  using  constitutive 
equations.  The  general  form  of  a  constitutive  equation  xmder  isothermal  conditions  for 
the  stress  tensor  as  a  function  of  the  strain  tensor  and  the  deformation-rate  tensor  may  be 
expressed  as  (Rohani  and  Thompson  1970) 

^rfij  Kn  ’  O  4.1 

where  is  the  strain  tensor  and  d„  is  the  deformation-rate  tensor.  The  response  function 
must  be  form  invariant  with  respect  to  rigid  motion  of  a  spatial  frame  of  reference  and 
invariant  with  respect  to  coordinate  transformations.  For  a  kinematically  linear  system, 
the  strain  rate  tensor  fejj  is  generally  assumed  to  be  the  same  as  the  deformation-rate 
tensor.  A  polynomial  of  two  symmetric  second  rank  tensor  variables  can  represent  in 

Equation  4.1  resulting  in  the  following  equation: 
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The  response  coefficients  tIq  through  qg  are  scalar- valued  functions  of  the  ten  joint 
invariants  of  Sy  and  fey  given  as 
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and  6y  is  the  Kroneker  delta.  The  response  coefficients  can  have  various  forms  and  must 
be  determined  from  experiments.  All  response  coefficients  are  not  required  for  all 
materials. 
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Based  on  the  information  presented  in  Chapter  3,  the  proposed  constitutive 
relationship  should  include  a  nonlinear  pressure-volume  relation,  nonlinear  shear  stress- 
strain  relation,  effect  of  pressure  on  the  shear  response,  failure,  and  fracture  as  well  as 
unloading  and  reloading  criteria.  A  mathematical  model  incorporating  these  response 
characteristics  will  include  expressions  of  the  form 

4.4 


or  the  hydrostatic  response,  and 
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for  the  deviatoric  response,  where  f,  and  fj  are  continuous  scalar  functions.  The 
volumetric  strain  and  octahedral  shear  strain  are  expressed  in  terms  of  the 
principal  invariants  of  strain  tensor  as 
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where  the  tensorial  shear  strains  (Sjj,  8,3,  and  823)  are  equal  to  one-half  of  the  engineering 
shear  strains  (Yu,  Y135  2nd  Y23)- 


Strain  rate  will  not  be  included  directly  in  the  development  of  the  current  model. 
It  is  clear  from  the  discussion  in  Chapter  3  that  brittle  geomaterials  are  affected  by 
loading  rate.  The  majority  of  available  data  however  are  concentrated  at  or  near  the 
unconfined  regime.  The  problems  of  interest  here  include  dynamic  loading  at  high 
pressures.  The  influence  of  loading  rate  as  confinement  increases  is  unclear.  The 
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available  information  on  the  effect  of  strain  rate  on  the  ultimate  strength  indicates  two 
approaches.  One  is  that  the  increase  in  dynamic  ultimate  strength  is  a  percentage  of  the 
static  ultimate  strength,  while  the  other  is  that  the  increase  in  dynamic  ultimate  strength  is 
essentially  a  shift  in  the  static  ultimate  strength  envelope  (that  is,  a  constant  increase  at  all 
pressure  levels).  An  indirect  approach  to  account  for  the  dynamic  effects  will  be  used 
here  by  increasing  the  static  ultimate  strength  envelope  by  a  percentage  without 
introducing  strain  rate  effects  explicitly. 


Eliminating  the  rate  effects  and  removing  the  higher  order  strain  terms  reduces 
Equation  4.2  to 

=  Tlo6..  +  T),e..  4.8 

The  response  coefficients  r^o  and  t],  must  now  be  selected  so  that  the  constitutive  equation 
yields  the  same  expression  for  mean  stress  and  octahedral  shear  stress  as  Equations  4.4 
and  4.5,  respectively.  Mean  stress  is  defined  as 
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Substituting  into  Equation  4.8  gives 
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As  stated  before, 
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where  Ij  is  the  second  invariant  of  the  deviatoric  strain  tensor  (Equation  4.7).  In  view  of 
Equations  4.8, 4.1 1  and  4.13, 
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Solving  Equation  4.14  for  qi  gives 
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4.15 


Substituting  into  Equation  4.10  and  solving  for  qo  gives 
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Substituting  Equations  4.15  and  4.16  into  Equation  4.8  and  rearranging  gives  the  general 
form  of  the  constitutive  equation  relating  the  stress  tensor  with  the  invariants  ”^oct>  Yoct» 
ande^  as 
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The  task  now  is  to  determine  appropriate  mathematical  relations  for  o„  and  /  Yo<.,. 


The  model  has  two  basic  premises:  (1)  the  material  response  can  be  decoupled 
into  hydrostatic  and  deviatoric  parts  and  (2)  the  deviatoric  response  can  be  further 
separated  into  cohesive  and  frictional  parts.  Assuming  the  response  to  be  decoupled 
implies  that  the  hydrostatic  (volumetric)  and  deviatoric  (shear)  parts  can  be  treated 
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separately.  In  general,  some  coupling  does  exist  since  the  response  of  a  material  to  pure 
shear  can  result  in  volumetric  deformation.  In  Equation  4. 1 7  the  term  to  the  left  of  the  + 
sign  represents  the  volumetric  part  and  the  term  to  the  right  of  the  +  sign  represents  the 
deviatoric  part.  A  hydrostatic  state  of  stress  will  result  in  volumetric  deformation  only 
and  a  state  of  stress  characterized  by  pure  shear  will  result  in  shear  deformation  only. 
Therefore,  each  part  can  be  fit  separately  to  the  volumetric  and  shear  response  of  the 
material. 


4.2.1  Volumetric  Response 


The  volumetric  response  defines  the  relationship  between  pressure  and  volumetric 
strain.  In  the  NIF  model,  this  response  is  separated  into  five  segments  as  a  function  of  the 
current  volumetric  strain  and  whether  loading,  unloading,  or  reloading  is  taking  place. 

The  volumetric  response  is  described  by  a  traision  portion  modeled  with  a  single  equation 
and  a  compression  portion  modeled  by  three  loading  equations  and  one  unload/reload 
equation.  Under  hydrostatic  tension,  the  response  of  the  material  has  the  same  initial 
modulus  Kg  as  the  compression  portion.  With  increasing  tension,  the  material  reaches  a 
maximum  level  of  stress  where  it  begins  to  soften  and  eventually  breakup.  The 
“loading”  response  in  tension  is  modeled  using  the  following  equation  developed  by  Elwi 
and  Murray  (1979) 


=  - 


1  + 


R  +  -  2 

K, 


X  -  {IR  -  \)X^  +  RX^ 


4.18 


where. 
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The  material  constants  are  8^,^,  8^  and  a^^and  are  shown  graphically  in  Figure  4.1. 
“Loading”  in  tension  occurs  if  (a)  the  volumetric  strain  is  negative  and  increasing  in 
magnitude,  or  (b)  an  unloading  from  a  compression  state  (positive  stress)  into  a  negative 
state  of  stress  and  the  volumetric  strain  is  decreasing.  No  data  is  available  on  which  to 
base  unloading  and  reloading  in  hydrostatic  tension.  “Unloading”  from  a  state  of  tension 
should  imply  an  increase  of  pressure.  Unloading/reloading  could  be  viewed  as  elastic  and 
follow  the  loading  response,  but  this  is  not  physically  realistic  once  softening  has  begun. 
The  unload/reload  could  also  be  “elastic”  by  following  a  straight  line  from  the  start  of  the 
unload/reload  to  zero  volumetric  strain.  This  approach  is  reasonable,  but  adds  another 
history  variable,  the  volumetric  strain  at  the  start  of  unloading,  that  must  be  tracked.  An 
intermediate  unloading  slope  could  be  used  but  it  is  not  known  what  the  response  should 
be  once  zero  pressure  is  reached.  The  last  option  is  to  assume  that  softening  does  not 
occur  during  hydrostatic  tension  such  that  the  pressure  remains  constant  once  is 
reached.  Unloading/reloading  is  then  treated  as  elastic  and  follows  the  loading  response. 
For  projectile  penetration  problems,  it  is  believed  that  the  details  of  the  unload/reload  in 
hydrostatic  tension  is  not  relevant  and  only  the  limiting  value  of  the  hydrostatic  tensile 
stress  (referred  to  as  the  tension  cutoff)  defined  by  is  important. 
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The  loading  and  unloading  portion  of  hydrostatic  compression  response  is  the 
same  as  that  described  in  Johnson,  et.  al.  (1994).  Initial  loading  in  compression  is  defined 
as 
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4.19 
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where,  8,,  is  the  volumetric  strain,  (o  is  the  pressure  at  {t^crmh  (8v)ot«a  » 

e;  =  8,  -  (eXci  and  K, ,  (8  a,  b,  c,  ,  and  (8^^*  are  material  constants. 
Unloading  in  the  compression  regime  occurs  if  (a)  the  volumetric  strain  is  decreasing  and 
(b)  it  is  less  than  the  previous  maximum  volumetric  strain.  Cmrently,  only  linear 

unloading/reloading  is  used.  If  the  unloading  begins  at  a  volumetric  strain  less  than 

» 

(8  the  unloading/reloading  follows  K^.  If  the  unloading  begins  at  a  volumetric 

strain  greater  than  (zXck  then  the  unloading/reloading  follows  If  the  unloading 
begins  at  a  volumetric  strain  between  {t^cnah  the  unloading/reloading  will 

follow  a  slope  given  by 
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If  the  unloading  extends  into  the  tension  (negative  pressure)  region  Equation  4.18  is  used 
to  define  the  response,  but  the  value  of  will  depend  on  where  the  unloading  originated. 
If  the  unloading  originated  at  a  strain  less  than  (8^)^^,  Equation  4.18  is  used  as  is.  If  the 
unloading  originates  at  a  strain  between  (e,,),^*  and  is  used  in  Equation  4.18 

instead  of  if  the  unloading  originates  at  a  strain  greater  than  (eXct^iock  is  used  in 
Equation  4. 1 8  instead  of  Each  section  of  the  compression  pressure-volumetric  strain 

response  and  the  corresponding  equation  are  shown  in  Figure  4.2. 
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4.2.2  Deviatoric  Response 


The  deviatoric  stress-strain  response  is  separated  into  friction  and  cohesion 
subparts.  The  octahedral  shear  stress  for  the  friction  subpart  increases  with  increasing 
strain  (strain  hardening)  and  is  a  function  of  the  superimposed  hydrostatic  state  of  stress. 
The  friction  subpart  is  described  by  an  equation  of  the  form  ~f  (Yoct»  O-  Th® 

following  hyperbolic  equation  is  used: 

^  ,  Iqc.  P  fa.  -  (^2  ^  ^^3 

f  (v,  -  v)  * 

where,  is  the  ultimate  strength  surface  containing  both  the  frictional  and  the  cohesive 
parts,  Xy^  is  the  cohesive  part  of  the  ultimate  strength  surface,  and  b^,  bj  and  define  the 
initial  modulus  of  the  x^^,  vs  y<,c/  response.  The  parameter  F  will  be  described  below.  The 
ultimate  strength  x^^  of  the  material  is  defined  using  a  two  part  curve  in  the  x^^,  vs 
space.  For  pressures  greater  than  zero  (compression),  the  curve  is  defined  as 

'^uit  =  ^  -  C  exp(BoJ  4.24 


For  pressures  less  than  zero  (tension),  the  curve  is  defined  as 


\it  =  (^  -  Q 


1.0  -  — 

o 

nicf 


4.25 


In  Equations  4.24  and  4.25,  A,  B  and  C  are  material  constants  that  define  the  nonlinear 
part  of  the  compression  side  of  the  ultimate  strength  surface.  Equation  4.25  is  the  linear 
part  of  the  tension  side  of  the  ultimate  strength  surface,  and  it  ensures  that  the  ultimate 
strength  of  the  material  at  the  maximum  hydrostatic  tensile  stress  is  zero.  Equations  4.24 
and  4.25  are  shown  graphically  in  Figure  4.3. 


For  the  cohesion  subpart,  the  stress  increases  with  increasing  strain  until  a  strain  is 
reached  where  the  maximum  cohesive  strength  of  the  material  is  exceeded.  This 
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Figure  4.3.  Ultimate  strength  envelope. 
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phenomenon  is  referred  to  as  cohesive  fracture.  Beyond  this  strain,  the  stress  will 
decrease  with  increasing  strain  (strain  softaiing).  The  cohesive  subpart  is  described  by  an 
equation  of  the  form  and  is  given  as 


T 


octc 


P  '^yc  Yc  ypct  ^ed 

^  V  Yc  ^  Yc  -  ^  y  yoct 


4.26 


where,  the  parameter  is  the  octahedral  shear  strain  at  which  strain  softening  begins, 
is  the  initial  slope  of  the  cohesion  subpart  of  the  deviatoric  stress-strain  response,  and  Rg^ 
is  a  parameter  that  defines  how  the  cohesive  fracture  will  progress,  i.e.,  for  very  brittle 
fracture  the  stress  will  decrease  rapidly  and  for  a  more  ductile  fracture  die  stress  will 
decrease  more  gradually.  The  parameter  y^  is  dependent  on  the  history  of  pressure  that 
has  been  applied  to  the  material.  The  maximum  cohesive  strength  is  dependent  on  the 
current  pressure  as  well  as  the  pressure  history  of  the  material.  As  was  discussed  in 
Chapter  3,  brittle  geomaterials  will  experience  internal  damage  during  the  application  of 
hydrostatic  pressure.  In  compression,  during  the  initial  phase  very  little  change  occurs. 
As  pressure  increases  the  level  of  “damage”  will  also  increase  as  bonds  between  the 
aggregates  are  broken.  Eventually  the  grains  will  rearrange  as  the  initial  cementation 
breaks  down.  Defining  the  boundaries  of  each  phase  is  based  on  the  different  regions  of 
the  pressure  vs  volumetric  strain  response  as  described  in  Figure  3.9.  The  onset  of  bond 
breakage  is  assumed  to  occur  when  (o  is  reached  using  Equation  4.19.  To  this  point 
in  compression,  is  constant  and  is 

Xyg=A-C  4.27 

which  is  the  value  of  at  equal  to  zero  from  Equation  4.24.  At  the  pressure  where 
the  cohesion  component  is  zero  (x^,,  is  zero),  the  response  is  only  fiictional.  This  point  is 
assumed  to  be  the  point  of  inflection  of  the  pressure  vs.  volumetric  strain  response  where 
it  begins  to  stiffen  rapidly.  The  value  of  (o is  defined  by  taking  the  second  derivative 
of  Equation  4.20,  setting  equal  to  zero  and  solving  for  e^.  This  value  of  (e^  =-b/3c)  is 
then  substituted  into  Equation  4.20  to  obtain 
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e  crush 


eL  + 

3c  27c2 


4.28 


Between  pressures  of  (o  and  (o  the  gradually  decreases  to  zero  using  the 
following  equation 


1.0  - 


m' crush 


^^rn^fric 


m  ’crush  / 
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The  ultimate  strength  envelope  of  the  cohesion  subpart  is  shown  graphically  in  Figure 
4.4. 


The  octahedral  shear  strain  at  which  softening  begins,  Yc»  should  be  reached  when 
the  total  calculated  octahedral  shear  stress  is  equal  to  the  ultimate  octahedral  shear 
stress  given  by  either  Equation  4.24  or  Equation  4.25.  Equations  4.23  and  4.26  are 
hyperbolic  equations  and  when  combined  provide  the  total  deviatoric  response  of  the 
material.  To  ensure  that  reaches  x„„  within  a  reasonable  value  of  (discussed  in 
Chapter  3),  a  factor  F  is  applied  in  Equations  4.23  and  4.26  and  is  described  below.  The 
value  for  Ye  where  equals  x„,,  can  be  determined  from  Equations  4.23  and  4.26  by  (1) 
substituting  Yc  for  Yoc/»  (2)  adding  the  two  equations,  (3)  setting  equal  to  1.0  since  we 
are  looking  for  the  value  of  Yc  where  strain  softening  just  begins,  (4)  setting  x„^,  equal  to 
x^ip  and  (5)  solving  for  Yc  to  give 


‘  {b2  -  Ki”.)"') 


4.30 


In  order  to  model  brittle  geomaterial  response  at  low  pressures  and  ductile  response  at 
high  pressures  with  a  transition  from  brittle  to  ductile  in  between,  the  following 
expressions  were  specified  for  F 


Equajtion  4.27  i 


Figure  4.4.  Maximum  strength  envelope  for  the  cohesion  subpart, 
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F  =  1.325 


F  = 


/ 

1.100  +  0.225 


1.000 


- 


m'crush 


m' crush  ) 


F  =  1.100 


for  o„  <  (o  4.31 

for  ^  o„  ^  (o„)^.,  4.32 

for  o„  >  (oj^.,  4.33 


The  following  equations  are  specified  fox 


=  1-0  forYoc/<Yc  4.34 

D  =  1.0  -  -1^ — for  Yc  ^  Yoct  ^  ^  Yc  4.35 
Yc  (^  -  1*0) 

i?ci  =  0.0  forYcc/>-^^Yc  4.36 

where  N  defines  how  far  Yoc/  is  firom  Yc  when  is  zero. 

Unloading  and  reloading  in  shear  will  be  determined  by  the  current  value  of  the 
octahedral  shear  strain.  If  is  less  than  the  largest  value  of  y^^,  obtained  thus  far  then 
the  material  is  assumed  to  be  either  unloading  or  reloading.  Load,  unload,  and  reload 
cycles  are  illustrated  in  Figure  4.5.  The  unload  response  initially  follows  a  slope  equal  to 
the  initial  total  modulus  of  the  vs  y^,  response  associated  with  the  pressure  at  the  time 

unloading  begins.  If  the  value  of  implied  by  the  total  modulus  and  the  current  value 
of  Yoc/  becomes  negative  (see  Figure  4.5),  then  the  response  is  determined  using 
Equations  4.23  through  4.33  with  =  1.0  so  that  no  softening  occurs.  The  value  for  Yo^ 
used  in  the  equations  is  the  difference  between  the  point  at  which  became  zero  and  the 

current  value  of  Yoc/-  Reloading  will  retrace  the  unloading  response  until  the  current  value 
of  Yoct  exceeds  the  maximum  value  of  Yoct  where  the  virgin  loading  of  the  material 


continues. 


edlAI  ‘ssejis  Jeaqs  leipeqepo 


Load/unload/reload  under  constant  pressure 
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The  model  that  has  been  developed  to  this  point  is  rate  independent.  From 
Figures  3.19  and  3.20  the  strength  of  conventional  concrete  increases  as  the  strain  rate 
increases.  The  ultimate  strength  in  shear  and  hydrostatic  tension  is  multiplied  by  a 
parameter  DYN  to  indirectly  account  for  rate  effect.  Currently,  DYN  is  a  constant  that 
increases  the  ultimate  strength  by  a  percentage. 


The  total  octahedral  shear  stress  is  obtained  by  adding  the  shear  stresses  from  the 
friction  (Equation  4.23)  and  cohesion  (Equation  4.26)  subparts  as 


oct 


“^octf  ^  “^octc 
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For  low  pressures  where  the  deviatoric  part  will  be  dominated  by  the  cohesive  subpart,  a 
brittle,  strain  softening  response  will  result.  As  pressure  increases,  the  frictional  subpart 
will  begin  to  dominate  and  the  response  will  become  more  ductile.  Figure  4.6  shows  the 
friction  and  cohesion  subparts  and  the  combined  response  for  a  constant  pressure  shear 
test.  The  stress  associated  with  the  friction  subpart  increases  monotonically  with  strain. 
The  stress  associated  with  the  cohesion  subpart  initially  increases  with  strain  but  begins 
to  decrease  beyond  Yo«  =  0.005.  The  combined  response  adds  both  subparts  to  result  in 
softening  behavior  and  a  residual  strength  (corresponding  to  the  frictional  subpart)  when 
the  cohesive  strength  is  completely  lost. 


The  NIF  model  requires  that  five  variables  be  stored  as  history  variables  for  each 
element  of  the  material.  These  variables  record  the  current  extremes  that  the  material  has 
experienced.  The  first  variable  is  the  maximum  volumetric  strain.  This  variable  is 
always  greater  than  or  equal  to  zero.  If  the  current  value  of  volumetric  strain  is  less  than 
this  value,  the  material  is  either  in  tension  or  in  an  unload/reload  cycle  with  respect  to 
pressure.  The  second  variable  is  the  maximum  octahedral  shear  strain.  This  value  will  be 
used  to  determine  if  the  material  is  in  an  unload/reload  cycle  with  respect  to  shear.  The 
third  variable  is  the  pressure  at  the  start  of  an  unload/reload  cycle  in  shear.  This  pressure 
will  be  used  to  determine  the  initial  total  modulus  of  the  vs  Yocr  response  when 
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unloading/reloading  in  shear.  The  fourth  variable  is  the  minimum  value  of  thus  far. 
The  fifth  variable  is  the  maximum  value  of  Yc  calculated  thus  far.  The  fourth  and  fifth 
variables  imply  that  healing  can  not  occur,  i.e.,  a  material  that  has  been  loading  to  the 
point  where  cohesion  has  been  lost  and  then  unloaded  does  not  regain  its  brittle 
characteristics  simply  because  it  has  been  unloaded. 

4.3  MODEL  FITTING 

The  NIF  model  must  now  be  fit  to  the  results  fi’om  laboratory  tests.  The  model 
has  nineteen  constants  that  are  defined  by  the  user:  eleven  constants  that  define  the 
pressure-volumetric  strain  response  (K^,  e^,  a„p  i£y)crush>  ^lock » ^nd  izJiock) 

and  eight  constants  that  define  the  x^,  vs  response  (bg,  bj,  b2.  A,  B,  C,  ,  and  JV).  The 
numerical  values  of  these  constants  must  be  determined  fi'om  laboratory  test  data  for  the 
material  of  interest.  To  illustrate  the  fitting  procedure,  the  model  will  be  fit  to  the 
conventional-strength  portland  cement  (CSPC)  concrete  whose  mix  proportions  are 
presented  in  Table  4. 1 . 

Constants  that  describe  the  tension  portion  of  the  pressure-volumetric  strain 
response  (K„  o^,  and  a„f)  can  be  taken  directly  fi-om  a  plot  similar  to  Figure  4.1. 

Very  little  data  is  available  for  the  hydrostatic  tension  response  of  brittle  geomaterials. 
Nichols  and  Ko  (1996)  present  data  for  a  plain  concrete  under  true  triaxial  tension 
(hydrostatic  tension).  Only  the  pre-peak  loading  portion  of  the  response  was  captured. 
The  average  pressure  at  failure  was  about  2.9  MPa  at  a  volumetric  strain  of  0.021  percent. 
They  note  that  the  stress  at  failure  is  the  same  as  that  fi-om  uniaxial  tension  experiments. 
Direct  pull  tension  experiments  on  CSPC  concrete  provided  an  average  strength  of  about 
3.6  MPa  at  an  axial  strain  of  about  0.012  percent.  Using  this  information  an  average 
value  of  3.4  MPa  will  be  used  for  The  value  for  K,  is  determined  fi-om  the 
compression  portion  of  the  response  that  is  described  below.  An  initial  bulk  modulus 
of  17,900  MPa  agrees  with  the  initial  slope  of  the  response.  Assuming  a  linear  response. 
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Table  4.1.  Ingredients  and  mixture  proportions  for  CSPC  concrete. 


Item 

Mixture  Proportions, 
Saturated  Surface-Dry 

Type  n  Portland  cement 

328.0  kg/m^ 

9.5-mm  limestone  coarse  aggregate 

1,034.1  kg/m^ 

Limestone  fine  aggregate 

806.3  kg/m^ 

Water 

186.9  kg/m^ 

Water  reducing  admixture 

1.3 1/m^ 

Air-detraining  agent 

0.33  kg/m^ 

w/c  =  0.57 
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a  value  for  of  0.019  percent  corresponds  to  the  above  values  for  and  K^.  Since  the 
response  becomes  nonlinear  as  the  peak  is  approached,  this  value  for  will  be  too  stiff, 
therefore,  a  value  of  0.030  percent  will  be  used.  Values  for  o^y-and  e^^in  Equation  4.18 
define  a  point  on  the  “descending”  portion  of  the  response.  Since  no  data  exists  on  which 
to  base  this  point,  a  value  of  aJ4  (0.9  MPa)  will  be  used  for  o„^and  4e^^  (0. 12  percent) 
will  be  used  for  Although  the  post-peak  softening  is  not  used,  it  is  recommended  to 
input  a^fSnA.  since  they  do  influence  the  pre-peak  portion.  The  response  from 
Equation  4.18  and  the  data  from  Nichols  and  Ko  are  shown  in  Figure  4.7. 

The  seven  constants  required  to  define  the  compression  part  of  the  pressure- 
volumetric  strain  response  are  (also  used  in  the  tension  part), 

(tj^ock  •  Th®  value  for  is  the  initial  slope  of  the  response  and  is  taken  to  be  17,900 
MPa.  The  value  for  (e  is  determined  as  the  point  on  the  pressure-volumetric  strain 
response  where  the  material  is  no  longer  linear  elastic  and  the  bonds  between  the 
aggregate  materials  begin  to  break  down.  For  the  CSPC  concrete,  a  value  of  0.229 
percent  will  be  used.  Values  for  and  (tJiocic  should  be  determined  next.  The 
recommended  pressure-volumetric  strain  response  shown  in  Figure  4.8  does  not  extend  to 
a  strain  where  the  modulus  is  nearly  constant.  Since  this  constant  modulus  should  occur 
when  all  of  the  air  has  been  compressed  out  of  the  material,  the  air  voids  content  should 
provide  some  assistance  in  determining  these  values.  The  locked  modulus  should  extend 
to  a  point  on  the  volumetric  strain  axis  that  is  close  to  the  air  voids  content.  From  the 
recommended  properties  in  Figure  4.6,  the  highest  unloading  modulus  is  74100  MPa. 
Extrapolating  this  modulus  from  a  volumetric  strain  equal  to  an  air  voids  content  of  6.4 
percent  shows  that  the  material  has  not  reached  the  locked  modulus.  The  pressure  at 
which  the  material  will  reach  the  locked  modulus  is  estimated  to  be  about  800  MPa.  This 
value  for  pressure  will  be  used  with  the  polynomial  portion  to  obtain  a  value  for  (tJiocic- 
Parameters  for  the  polynomial  part  of  the  response  between  (eJcnuA  i^J\ock 
determined  with  the  aid  of  curve  fitting  software.  Values  for  a,  b,  and  c  are  6244.15  MPa, 
-47387.5  MPa,  and  1212580.0  MPa,  respectively.  The  value  for  is  obtained  by 
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determining  the  volumetric  strain  that  coincides  with  a  pressure  of  800  MPa,  and  was 
calculated  to  be  7.98  percent.  The  model  fit  and  recommended  response  are  shown  in 
Figure  4.9. 

The  recommended  properties  for  the  CSPC  concrete  do  not  include  experiments 
conducted  imder  constant  pressure,  which  ideally  is  needed  to  fit  the  NIF  model.  The  set 
of  recommended  properties  does  include  the  results  from  constant  radial  stress 
experiments,  which  can  be  used  if  the  ultimate  strength  envelope  is  assumed  to  be  path 
independent,  i.e.,  the  same  envelope  is  obtained  regardless  of  the  loading  path  during  the 
experiment.  Constants  that  describe  the  deviatoric  part  of  the  model  (Jb^,  bj,  b2,  Pc' 

C,  and  N)  are  obtained  from  the  results  of  shear  experiments.  Values  for  A,  B,  and  C  are 
obtained  by  fitting  Equation  4.24  to  the  peak  strengths  from  the  shear  experiments  plotted 
as  vs  o„.  Using  the  curve  fitting  software,  values  for  ^4,  B,  and  C  were  determined  to 
be  165.197  MPa,  -0.0044034265  MPa'*,  and  156.913  MPa,  respectively.  A  comparison 
of  the  model  ultimate  strengths  (total,  cohesion,  and  fiiction)  and  the  recommended 
response  is  shown  in  Figure  4.10.  The  total  envelope  from  the  model  agrees  very  well 
with  the  recommended  surface.  The  cohesion  part  follows  the  total  envelope  in  the 
negative  pressure  region,  is  equal  to  the  octahedral  shear  stress  intercept  from  zero 
pressure  until  is  reached,  and  then  decreases  to  zero  at  (o  The  friction  part  is 

the  difference  between  the  total  envelope  and  the  cohesion  part. 

Values  for  bg,  bj,  and  are  determined  by  fitting  the  shear  modulus  part  of 
Equation  4.23  to  the  initial  shear  modulus  from  the  results  of  shear  experiments 
conducted  at  several  pressure  levels.  The  data  provided  with  the  recommended  properties 
can  be  used  by  noting  that  the  shear  modulus  can  be  obtained  from  the  plots  of  principal 
stress  difference  vs  principal  strain  difference.  Values  selected  for  bg,  bj,  and  b2  are 
2530.21, 0.192653,  and  2758.62  MPa,  respectively.  A  comparison  of  the  model  shear 
modulus  and  the  data  is  presented  in  Figure  4. 1 1 .  The  shear  modulus  for  the  cohesion 
part  Pc  cannot  be  obtained  from  the  data  provided  in  the  recommended  properties  since 
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fit  and  recommended  ultimate  strength  envelope. 


15,000 


100 


edl/\l  'sninpoiAi  jeaqs 


Figure  4.1 1.  NIF  model  shear  modulus  and  values  obtained  from  the  recommended  CSPC  concrete  properties. 
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the  friction  part  is  included  in  all  experiments  for  which  recommended  properties  were 
provided.  Using  the  shear  modulus  provided  for  the  experiment  at  zero  confimng 
pressure  (14000  MPa)  should  approximate  this  value.  The  last  parameter  required  for  the 
shear  part  of  the  model  is  N.  The  value  of  N  should  be  based  on  the  softening  portion  of  a 
pure  shear  experiment  conducted  at  zero  pressure,  for  which  little  or  no  data  exists.  The 
recommended  properties  do  not  contain  such  data.  An  extrapolation  for  the  softening 
portion  of  a  dynamic  imconfined  compression  experiment  provided  with  the  properties 
indicates  that  the  strength  is  zero  at  about  an  of  4.  For  lack  of  a  better  value,  4  will  be 
used  unless  results  from  exercising  the  model  indicate  a  change  is  needed.  The  nineteen 
constants  for  the  NIF  model  fit  to  the  CSPC  concrete  are  summarized  in  Table  4.2. 

A  second  concrete  that  was  used  in  several  projectile  penetration  experiments  was 
also  fit  to  the  NIF  model.  The  mix  proportions  for  the  WES5000  concrete  are  presented 
in  Table  4.3.  The  procedure  used  to  fit  the  model  to  this  concrete  are  the  same  as  those 
used  to  fit  the  model  to  the  CSPC  concrete.  The  nineteen  constants  for  the  NIF  model  fit 
to  the  WES5000  concrete  are  summarized  in  Table  4.4. 

4.4  MODEL  DRIVER 

The  NIF  model  has  now  been  “fit”  to  a  set  of  recommended  properties  by 
determining  the  numerical  values  for  each  of  the  model  constants.  The  next  step  in 
determining  the  adequacy  of  the  model  and  the  quality  of  the  fit  is  to  use  a  computer  code 
containing  the  model  (referred  to  as  a  model  driver)  to  drive,  or  exercise,  the  model  along 
loading  paths  to  which  the  model  has  been  fit,  and  paths  which  can  demonstrate  the  level 
of  robustness  of  the  model.  The  driver  can  also  be  used  to  make  adjustments  to  the 
model  parameters  to  provide  better  agreement  with  the  recommended  properties.  If  a 
particular  loading  path  has  been  found  to  be  prominent  in  the  problem  for  which  the 
model  will  be  used,  adjustments  can  be  made  to  appropriate  model  parameters  to  better 
replicate  that  path.  The  response  along  other  loading  paths  may  be  sacrificed,  if 
necessary,  in  doing  this  and  must  be  checked  for  adequacy. 


Table  4.2.  Summary  of  NIF  model  constants  for  the  CSPC  concrete. 


NIF  Model  Constant 

Constant  Value 

17,900  MPa 

Sv. 

0.0003 

V 

0.0012 

me 

3.4  MPa 

V 

0.9  MPa 

^^ervsh 

0.00229 

a 

6,244.15  MPa 

b 

-47,387.5  MPa 

c 

1,212,580.0  MPa 

^lock 

74,100  MPa 

0.0798 

bo 

2,530.21 

b, 

0.192653 

b: 

2,758.62  MPa 

A 

165.197  MPa 

B 


-0.0044034265  MPa'’ 
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Table  4.3.  Ingredients  and  mixture  proportions  for  WES5000  concrete. 


Item 

Mixture  Proportions, 
Saturated  Surface-Dry 

Type  I  Portland  cement 

264.0  kg/m^ 

Fly  ash 

55.8  kg/m^ 

9.5-mm  local  unprocessed  chert  coarse  aggregate 

1,037.6  kg/m^ 

Local  unprocessed  chert  fine  aggregate 

840.7  kg/m^ 

Water 

145.9  kg/m^ 

Water  reducing  admixture  “300N” 

0.65 1/m^ 

High-range  water  reducing  admixture  “Rheobuild  716" 

1.6 1/m^ 

w/c  =  0.46 


Table  4.4.  Summary  of  NIF  model  constants  for  the  WES5000  concrete. 
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The  driver  is  a  computer  program  that  runs  on  a  PC  that  can  exercise  the  model 
along  any  prescribed  stress  or  strain  load/unload/reload  path.  An  undocumented  program 
obtained  from  Akers  (1992)  was  modified  for  this  purpose.  The  model  has  been 
developed  to  accept  strains  as  input  in  order  to  obtain  the  corresponding  stresses.  The 
driver  has  been  set  up  so  that  either  stress  or  strain  paths  can  be  prescribed. 

The  total  strain  based  on  the  original  configuration  of  the  material  will  be  supplied 
to  the  model.  The  sequence  of  operation  of  the  model  is  to  first  calculate  the  current 
pressure  based  on  the  current  volumetric  strain.  In  the  driver,  the  volumetric  strain  is 
calculated  using  Equation  4.6.  The  criteria  described  above  in  Section  4.3  are  followed  to 
determine  whether  the  material  is  loading,  unloading  or  reloading.  Once  the  pressure  has 
been  determined,  the  deviatoric  part  is  calculated  following  the  criteria  for  loading, 
unloading,  and  reloading  described  above.  Equation  4. 1 7  is  then  used  to  determine  the 
individual  stress  components.  The  coding  of  the  subroutine  containing  the  NIF  model 
that  is  used  in  EPIC  (provided  in  the  Appendix)  is  essentially  the  same  as  the  coding  used 
in  the  driver. 

A  comparison  of  the  model  driver  and  the  recommended  pressure- volumetric 
strain  response  for  the  CSPC  concrete  is  shown  in  Figure  4.12.  Results  from  the  model 
driver  compare  well  with  the  recommended  response  for  initial  loading.  The  unloadings 
from  the  model  driver  do  not  agree  as  closely  since  the  model  currently  has  no 
mechanism  to  allow  for  the  hysteresis  that  occurs  as  the  pressure  approaches  zero.  A 
comparison  of  the  model  driver  and  the  recommended  response  for  constant  confining 
pressure  triaxial  shear  (CTC)  experiments  is  shown  in  Figure  4.13.a  as  principal  stress 
difference  (o^j^)  vs  principal  strain  difference  (e^^)  and  Figure  4.13.b  as  principal  stress 
difference  vs  axial  strain  during  the  shear  phase  where 


^diff  ~ 


4.38 
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Figure  4.12.  Comparison  of  NIF  model  driver  and  recommended  response  for  hydrostatic  compression, 
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and 

W  '  fh  4.39 

The  model  driver  agrees  reasonably  well  with  the  recommended  responses.  Figure  4. 1 3 
also  illustrates  the  capacity  of  the  model  to  simulate  the  transition  from  brittle  to  ductile 
response. 

A  comparison  of  the  model  driver  response  and  the  experiment  results  from 
Figure  3.12  are  shown  in  Figure  4.14  for  a  uniaxial  state  of  stress.  The  model  agrees  well 
with  the  axial  stress  versus  lateral  and  axial  strain  to  peak  stress.  The  model  also  shows 
hi^er  values  of  axial  strain  than  radial  strain  similar  to  the  experiment  results.  The 
model  response  then  softens,  but  the  softening  portion  from  the  experiment  was  not 
recorded.  The  axial  stress  versus  volumetric  strain  from  the  model  is  stiffer  than  the 
experiment  results  and  does  not  show  a  tendency  to  dilate  near  the  peak.  One  explanation 
for  this  response  is  that  the  experiment  is  less  “controlled”  than  the  driver  in  that  the 
driver  adjusts  the  lateral  strain  in  an  attempt  to  maintain  zero  radial  stress. 

An  additional  loading  path  provided  with  the  recommended  responses  is  a 
uniaxial  strain  loading  in  compression.  A  comparison  between  the  model  driver  and  tiie 
recommended  response  is  shown  in  Figure  4.15.  The  axial  stress  vs  axial  strain  responses 
agree  well,  but  the  principal  stress  difference  vs  mean  normal  stress  (pressure)  response 
from  the  model  driver  is  slightly  below  the  recommended  response.  Unloading  imder 
uniaxial  strain  conditions  results  in  unloading  in  both  pressure  and  shear.  This  results  in 
some  hysteresis  within  the  model  since  the  shear  modulus  is  changing  as  the  pressure 
changes,  particularly  once  the  principal  stress  difference  reaches  zero.  The  model 
response  shows  continued  unloading  resulting  in  the  “J-shape”  of  the  stress  path  in  the 
negative  principal  stress  difference  regime  (obtainable  since  the  specimen  is  assumed  to 
be  axisymmetric).  Data  to  support  this  response  is  not  available  for  CSPC  concrete,  but 
similar  responses  have  been  observed  for  clayey  sands  (Figure  4.16). 
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Figure  4.14.  Comparison  of  NIF  model  driver  and  experiment  results  for  an  unconfined  compression  experiment. 
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Figure  4.15.  Comparison  of  NIF  model  driver  and  recommended  responses  for  uniaxial 
strain  loading. 
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Although  not  used  to  fit  the  model,  constant  pressure  (TC)  experiments  were 
conducted  on  specimens  of  the  WES5000  concrete.  Results  from  experiments  and  the 
model  driver  at  constant  pressures  of  25  MPa  after  pre-stresses  of  0, 100,  and  400  MPa 
are  shown  in  Figure  4. 1 7.  The  experiment  results  show  that  a  pressure  of  25  MPa  is 
below  the  crush  pressure  for  these  specimens,  100  MPa  is  close  to  the  crush  pressure,  and 
400  MPa  is  well  beyond  the  crush  pressure  and  may  be  close  to  or  beyond  the  fiiction 
pressure.  Results  from  pre-stress  experiments  of  0  and  100  MPa  are  similar.  But  the  pre¬ 
stress  of 400  MPa  resulted  in  a  softer  shear  response  and  slightly  lower  strength.  The 
model  driver  indicates  the  pre-stress  of  100  MPa  is  beyond  the  crush  pressure  and  400 
MPa  is  beyond  the  friction  pressure.  The  model  driver  therefore  results  in  a  softer  shear 
response  as  the  pre-stress  increases,  and  a  slightly  higher  ultimate  strength  since  the  value 
of  Yc  increases  as  pressure.  Unlike  the  experiment  results,  the  model  driver  does  not 
indicate  dilation  since  the  test  boundary  conditions  imply  zero  change  in  volumetric  strain 
for  constant  pressure.  Results  from  experiments  and  the  model  driver  at  a  constant 
pressure  of  100  MPa  after  pre-stresses  of  0  and  400  MPa  are  shown  in  Figure  4.18.  The 
experiment  results  show  a  slightly  stiffer  response  for  the  pre-stressed  material.  The 
pressure- volumetric  strain  response  for  the  pre-stressed  material  is  stiffer  than  the 
response  for  the  non-prestressed  material,  which  may  explain  the  difference.  The  model 
driver  shows  a  softer  response  for  the  pre-stressed  material,  which  results  from  the  loss  of 
cohesion  and  increase  in  Yc  as  the  pressure  increased,  and  the  assumption  that  unloading 
does  not  “heal”  the  material. 
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Volumetric  Strain,  percent  Axial  Strain  During  Shear,  percent 

b.  NIF  model  driver 

Figure  4.17.  Comparison  of  NIF  model  driver  and  experiment  results  for  TC  experiments  at  25  MPa  after  pre-stresses  of  0, 100, 
and  400  MPa. 
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18.  Comparison  of  NIF  model  driver  and  experiment  results  for  TC  experiments  at  100  MPa  after  pre-stresses  of  0 
and  400  MPa. 
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CHAPTER  5 

ftnttf-ft  .fmf.nt  code  implementatton  and  evaluation 

5.1  FINITE-ELEMENT  CODE  IMPLEMENTATION 

The  constitutive  model  developed  in  Chapter  4  must  now  be  implemented  into  a 
numerical  computer  code  in  order  to  be  a  useful  tool.  The  computer  code  selected  for  this 
purpose  is  the  EPIC  code  (Johnson,  et.  al.  1994).  This  process  involves  flie 
implementation  of  algorithms  to  calculate  the  correct  strain  input  to  the  model  and  to 
convert  the  stress  output  by  the  model  into  a  compatible  form  accepted  by  the  code.  The 
NIF  model  has  been  developed  based  on  the  assumption  that  flie  strain  input  to  the  model 
is  the  total  strain  related  to  the  original  configuration  (total  change  in  lengtii/original 
length)  and  the  output  stress  is  the  Cauchy  stress  aligned  with  the  original  configuration. 
The  output  stress  will  need  to  be  rotated  to  the  current  configuration  of  the  element  to  be 
compatible  with  the  kinematics  of  the  code.  Other  issues  that  will  be  addressed  are  the 
use  of  artificial  viscosity  and  criteria  for  elimination  of  highly  distorted  elements  firom  a 
simulation. 

The  EPIC  code  is  capable  of  solving  two-  and  three-dimensional  projectile  impact 
and  penetration  problems.  The  code  has  several  element  options  and  a  material  model 
library  for  several  materials  typically  encountered  in  these  types  of  problems.  All  of  the 
available  models  for  geomaterials  are  based  on  plasticity  concepts.  The  code  can  account 
for  the  sliding  interface  between  projectile  and  target  and  can  erode  highly  deformed 
elements.  The  computational  flow  of  the  code  is  shown  schematically  in  Figure  5.1 
(Johnson,  et.  al.  1978).  The  geometry  of  the  problem  is  first  represented  with  elements. 
The  mass  of  the  material  contained  in  each  element  is  then  lumped  at  the  nodes  and 


Initial  Conditions  Integration  Loop 


Figure  5.1.  Computational  technique  used  in  EPIC  (Johnson,  et  al.  1978). 
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velocities  are  assigned.  The  velocities  are  used  to  obtain  the  displacements  and  the 
strains  and  strain  rates  are  calculated.  The  stresses  within  the  element  are  then  calculated 
using  the  mathematical  material  model.  The  stresses  are  converted  into  concentrated 
forces  at  the  nodes  and  the  nodal  accelerations  are  calculated  using  Newton’s  equations  of 
motion.  The  equations  of  motion  are  evaluated  to  determine  the  updated  velocities  at  the 
nodes  and  the  process  is  repeated  until  a  termination  criterion  (maximum  time)  is  met.  In 
this  scheme,  the  NIF  model  will  be  used  in  the  computation  of  the  stresses  within  the 
element  consisting  of  a  brittle  geomaterial.  Models  used  for  computing  stresses  within 
the  projectile  will  not  be  altered. 

5.1.1  Calculation  of  Strain  and  Stress 

The  EPIC  code  is  based  on  an  explicit  Lagrangian  finite-element  formulation 
where  the  equations  of  motion  are  integrated  directly.  This  eliminates  the  need  for  the 
traditional  stiffness  matrix  approach.  The  general  form  of  the  code  uses  an  updated 
Lagrangian- Jaumann  kinematic  formulation  to  account  for  the  finite  strains.  This 
formulation  uses  the  Cauchy  stress  (o)  and  the  rate  of  deformation  (/)),  respectively,  as 
the  stress  and  strain  measure,  and  the  Jaumann  stress  rate  tensor  as  the  objective  stress 
rate.  The  frame  of  reference  is  the  current  configuration  of  the  body. 

The  total  Lagrangian  kinematic  formulation  (Adley,  et.  al.  1996)  that  will  be  used 
to  implement  the  NIF  model  into  EPIC  will  use  the  Green-Lagrange  strain  as  input  to  the 
model.  The  non-conjugate  Cauchy  stress  associated  with  the  original  orientation  (o^)  will 
be  ouq)ut  by  the  model  (Bazant  1996  and  Bazant,  et.  al.  1996). 

The  Green-Lagrange  strain  required  as  input  to  the  NIF  model  will  be  computed 
as  follows  using  the  rate  of  deformation  tensor  (/))  and  the  spin  tensor  (fV)  that  are 
calculated  by  the  code.  The  spatial  velocity  gradient  (X)  is  related  to  the  deformation 
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gradient  (F)  by 


F  =  L  F 


5.1 


where 


L  =  D  +  W 


5.2 


Assuming  that  L  is  constant  during  the  current  time  step.  Equation  5. 1  has  as  a  solution  at 
the  current  time  (t) 

m  =  exp[i(<  -  g]  5.3 

where  F„  is  the  value  of  F  at  the  end  of  the  preceding  time  step  (t„).  The  exponential 
function  in  Equation  5.3  can  be  approximated  by  the  Taylor  series  expansion  as 

exp[£(r  -  O]  =  I  +  L(t  -  t„)  +  -  tf  +  ...  5.4 


where  I  is  the  identity  matrix  and  higher  order  terms  are  not  considered  since  the  time 
increment  is  very  small  for  the  problems  of  interest.  The  increment  of  F  during  the 
current  time  step  (AF)  is  then  calculated  as 


AF  =  F(t)  -  Fit„)  = 


LAt  +  -L^At^ 
2 


5.5 


where  At  =  t- 1„.  The  Green-Lagrange  strain  (e)  at  the  current  time  can  then  be  calculated 
as 


e  =  }-(f^F  -  5.6 

2 

with  F  =  F„  +  AF  where  F„  is  the  deformation  gradient  at  the  previous  time.  Note  that 
the  off-diagonal  terms  in  e  are  the  tensorial  values  of  the  shear  strain,  i.e.,  e  j2  =  Y12  /  2. 
These  strains  will  be  used  in  the  NIF  model  to  calculate  the  Cauchy  stresses  associated 
with  the  original  configuration  (0''). 
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The  stresses  calculated  by  the  NIF  model  must  be  rotated  back  to  the  current 
configuration  before  continuing  within  the  EPIC  code.  The  transformation  used  is 

o  =  Ro^R^  5.7 


The  rotation  matrix  (R)  can  be  determined  from  information  contained  within  EPIC  using 
an  equation  based  on  the  relationship  between  F  and  the  left-hand  stretch  tensor  (V)(F= 
VR)  developed  by  Anderson,  et.  al.  (1994)  as 


/  +  A/ 


(  1  1 
/  -  -AtW 

-1 

<1 

1 

+ 

2  J 

1  2  ; 

5.8 


The  spin  tensor  (W)  can  be  obtained  in  EPIC  fi"om  the  nodal  velocity  (v)  as 


W.J  =  - 
2 


( 


dv.  dv.  ^ 
dyj  dyj  j 


5.9 


The  updated  rotation  tensor  from  Equation  5.8  is  used  in  Equation  5.7  to  calculate  the 
Cauchy  stress  related  to  the  current  configuration  of  the  element. 


5.1.2  NIF  Model  Subroutine 


The  subroutine  (NIF)  used  to  implement  the  NIF  model  into  the  EPIC  code  is 
presented  in  the  Appendix.  The  first  part  of  the  subroutine  consists  of  the  DIMENSION 
and  PARAMETER  statements  required  by  the  subroutine.  Constants  that  will  be  used  in 
the  subroutine  and  the  NIF  model  constants  are  then  assigned  to  the  code  variables.  All 
of  the  code  variable  names  assigned  to  the  material  constants  begin  with  an  “A”  only  as  a 
matter  of  convenience.  The  next  section  of  the  subroutine  computes  the  total  stresses  and 
stress  increments  from  the  previous  time  for  use  in  an  energy  calculation  later  in  the 
subroutine.  The  next  section  of  the  subroutine  uses  the  procedures  outlined  in  Section 
5.1.1  to  calculate  the  Green-Lagrange  strains  and  the  rotation  matrix. 
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The  volumetric  strain  input  to  the  model  for  calculation  of  the  pressure  is  the 
volumetric  strain  calculated  within  the  EPIC  code  based  on  the  deformation  of  the 
element.  An  additional  variable  computed  as  part  of  the  model  is  EBAR.  This  variable 
will  be  used  to  determine  when  an  element  has  deformed  to  an  extent  that  it  no  longer  can 
provide  meaningful  resistance  to  the  applied  deformations,  and  may  be  controlling  the 
computational  speed  of  the  simulation  due  to  its  small  aspect  ratio.  The  variable  is 
calculated  in  shear  as  the  ratio  of  the  current  value  of  Yoc<  /  Yc  •  Since  failure  can  also 
occur  in  hydrostatic  tension,  a  similar  quantity  is  calculated  as  .  The  maximum  of 
the  two  quantities  is  assigned  to  EBAR.  EBAR  in  essence  tracks  how  far  beyond  the 
critical  strain  at  the  start  of  softening  the  material  has  strained.  EBAR  is  checked  against 
a  user  input  (EBAR^riO  ^  determine  if  the  element  is  to  be  removed  from  the  calculation 
by  setting  all  stresses  to  zero,  but  retaining  the  mass  of  the  element  at  the  nodes.  An 
advantage  to  this  approach  is  that  EBAR  becomes  dependant  on  the  current  state  of  stress 
within  the  element  and  is  not  simply  a  scalar  value  that  checks  the  state  of  plastic  strain 
within  the  element,  which  is  commonly  used.  At  high  pressmes,  the  use  of  plastic  strain 
as  the  check  may  prematurely  eliminate  elements  that  are  still  providing  a  sigmficant 
resistance  to  deformation.  The  above  calculation  of  EBAR  will  result  in  a  criterion  at 
high  pressures  that  is  larger  than  the  criterion  at  low  pressures  which  makes  more  sense 
from  a  physical  standpoint. 


An  artificial  viscosity  {Q)  is  also  used  with  the  model.  The  artificial  viscosity  is 
applied  to  the  normal  stresses  to  damp  out  localized  oscillations  of  the  concentrated 
masses  that  might  otherwise  occur  during  wave  propagation  problems.  The  approach 
used  is  the  original  method  contained  in  the  EPIC  code  (Johnson,  et.  al.  1997).  The 
equation  for  Q  contains  linear  and  quadratic  components  and  is  expressed  as 


Q  = 


+  CgphHef 


5.10 


where  and  Cq  are  the  linear  and  quadratic  dimensionless  coefficients,  respectively,  p  is 
die  density,  is  the  soxmd  velocity  of  the  material,  and  h  is  the  minimum  altitude  of  the 
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element.  Equation  5.10  is  only  applied  if  the  volumetric  strain  rate  (ej  is  compressive. 
The  value  of  0  is  a  pressure  that  is  added  to  the  total  normal  stresses  output  from  the  NIF 
model  when  the  element  is  being  compressed.  When  the  element  is  expanding,  no 
artificial  viscosity  is  applied.  Values  of  1.75  for  Q  and  14.0  for  Cq  give  smooth  results 
for  one-dimensional  wave  propagation  calculations  using  the  model. 

5.2  MODEL  IMPLEMENTATION  EVALUATION 

The  NIF  model  as  implemented  into  the  EPIC  code  must  be  evaluated  to  ensure 
that  the  implementation  5delds  the  desired  results.  This  evaluation  will  consist  of 
exercising  the  model  for  a  single  element  along  prescribed  boundary  conditions  to 
simulate  a  hydrostatic  compression  loading,  an  unconfined  compression  loading,  and  a 
uniaxial  strain  compression  loading.  Unload/reload  cycles  will  also  be  included  with 
some  cycles  proceeding  into  extension. 

The  element  configuration  and  loading  history  for  hydrostatic  compression  are 
shown  in  Figure  5.2.  The  simulation  is  for  a  single  axisymmetric  quad  element  with  the 
integration  point  at  the  center  of  the  element.  Pressure  applied  to  the  three  outer  sides  of 
the  element  varies  with  time  to  provide  several  unload/reload  cycles  to  a  maximum 
pressure  of  600  MPa  at  a  time  of  0.00245  s.  Results  from  the  simulation  are  compared  to 
the  results  from  a  similar  loading  using  the  NIF  model  driver  in  Figure  5.3.  Both  the 
EPIC  simulation  and  the  NIF  model  driver  provide  the  same  responses.  Results  from  the 
EPIC  simulation  also  show  some  of  the  unload  cycles  extending  into  the  hydrostatic 
tension  region  where  the  pressure  is  limited  to 

The  element  configuration  and  loading  history  for  unconfined  compression  are 
shown  in  Figure  5.4.  The  element  configuration  is  the  same  as  that  for  hydrostatic 
compression.  The  loading  is  applied  by  giving  the  nodes  a  compressive  velocity  equal  to 
101.6  mm/s  over  the  duration  of  the  simulation.  Results  from  the  EPIC  simulation  are 
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Figure  5.2.  Single  element  configuration  and  loading  history  used  in  EPIC 
to  check  the  NIF  model  for  hydrostatic  compression. 
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Figure  5.3.  Comparison  of  the  hydrostatic  compression  from  EPIC  and  the  NIF  model 
driver. 
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Figure  5.4.  Single  element  configuration  and  loading  used  in  EPIC 
to  check  the  NIF  model  for  unconfined  compression. 
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compared  to  the  results  using  the  NIF  model  driver  in  Figure  5,5.  Except  for  the  point 
where  the  axial  stress  has  reduced  to  zero  during  the  softening  phase,  the  results  are  very 
similar.  The  discrepancy  at  the  end  of  the  softening  is  believed  to  come  fi'om  two 
sources.  First  is  a  small  contribution  by  the  artificial  viscosity  used  in  the  EPIC 
simulation.  Second  is  that  the  element  is  fi-ee  to  deform  solely  based  on  the  resistance 
supplied  by  the  model,  but  the  applied  velocity  is  constant  in  the  vertical  direction.  This 
may  act  to  introduce  some  “confinement”  to  the  element. 

The  element  configuration  and  loading  history  for  uniaxial  strain  are  shown  in 
Figure  5.6.  The  element  configuration  is  again  a  single  axisymmetric  quad,  but  the 
bottom  nodes  are  constrained  from  moving  and  the  top  nodes  are  constrained  in  the 
lateral  direction.  The  loading  is  applied  by  giving  the  top  nodes  a  varying  compressive 
velocity  to  provide  unload/reload  cycles  to  a  maximum  velocity  equal  to  30.4  m/s. 
Results  fi-om  the  EPIC  simulation  are  compared  to  the  results  using  the  NIF  model  driver 
in  Figure  5.7.  Again,  results  fi-om  the  EPIC  simulation  agree  well  with  the  results  fi*om 
the  NIF  model  driver. 

Based  on  the  above  simulations,  the  NIF  model  appears  to  be  properly 
implemented  into  the  EPIC  code  and  fimctioning  as  expected.  In  the  following  chapter, 
simulations  will  be  performed  to  evaluate  using  the  NIF  model  in  EPIC  to  simulate 
projectile  penetration  and  perforation  experiments. 
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Figure  5.5.  Comparison  of  unconfined  compression  from  EPIC  and  the  NIF  model  driver. 
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Figure  5.6.  Single  element  configuration  and  loading  used  in  EPIC 
to  check  the  NIF  model  for  uniaxial  strain. 
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Figure  5.7.  Comparison  of  uniaxial  strain  loading  from  EPIC  and  the  NIF  model  driver. 
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CHAPTER  6 

PENETRATION  EXPERIMENTS  AND  NUMERICAL  SIMULATIONS 

6.1  PENETRATION  EXPERIMENTS 

Performance  of  the  NIF  model  in  EPIC  will  be  evaluated  by  comparing  code 
simulations  to  the  results  from  projectile  penetration  and  perforation  experiments.  The 
experiments  involved  the  launch  of  steel  projectiles  into  concrete  targets  in  the  WES  83- 
mm  ballistic  range.  Penetration  experiments  were  conducted  using  a  robust,  thick-walled 
armor  piercing  (AP)  projectile  launched  into  semi-infinite  (target  thickness  much  greater 
than  the  penetration  depth)  CSPC  concrete  targets  at  impact  velocities  ranging  from  277 
to  800  m/s.  Perforation  experiments  were  conducted  using  a  robust,  semi-armor  piercing 
(SAP)  projectile  launched  into  WES5000  concrete  targets  witii  thicknesses  ranging  from 
127  to  284  mm.  In  both  series  of  experiments  the  projectile  sustained  minimal  damage. 
The  following  sections  describe  the  WES  83-mm  ballistic  range  and  the  results  from  the 
experiments. 

6.1.1  WES  83-mm  Ballistic  Range 

The  WES  83-mm  Ballistic  Range  (Frew,  et.al.  1993)  consists  of  an  83-mm  solid- 
propellent  launcher,  a  mount  to  support  and  align  the  launcher,  a  blast  tank,  a  sabot 
separator  system,  a  drift  tube  assembly,  and  a  target  room  as  shown  in  Figure  6. 1 .  The 
83-mm  "gun"  is  capable  of  launching  projectiles  with  masses  up  to  2.8  kg  at  velocities  in 
excess  of  2  km/s  and  projectiles  with  masses  of  12  kg  or  more  at  velocities  of  1  km/s.  At 
the  downstream  end  of  the  launch  tube  a  vented  section  extends  into  the  blast  tank.  The 
vents  (large  ports  cut  into  the  tube  wall)  allow  the  accelerating  gases  to  expand  laterally 
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Figure  6.1.  WES  projectile  penetration  research  facility. 
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into  the  blast  tank.  The  effect  of  the  vent  section  is  to  terminate  projectile  acceleration 
and  muzzle  the  sound.  The  blast  tank  contains  two  central  baffles  used  to  break  up  large 
gas  flows  and  is  equipped  with  a  ventilation  system  to  remove  the  explosive  gas  after  the 
test.  The  projectile  travels  down  the  launch  tube  in  a  hard  plastic  sabot.  The  sabot  can 
either  be  stripped  using  a  sabot  separator  system,  stripped  aerodynamically  or  left  on  the 
projectile.  The  sabot  separator  system  consists  of  a  gasdynamic  tube  and  impact  tank. 
The  gasdynamic  tube  is  an  extension  of  the  launch  tube  with  provisions  to  seal  each  end 
of  the  tube.  The  gasdynamic  tube  can  be  operated  at  atmospheric  pressure  or  pressurized 
at  levels  up  to  2.0  atm.  A  sabot  impact  tank  is  used  to  intercept  and  pulverize  oncoming 
sabots  while  allowing  projectiles  to  pass  through  the  drift  tube  and  into  the  target  room 
unimpeded.  If  stripped  aerodynamically,  the  multi-piece  sabot  folds  away  from  the 
projectile  as  it  is  traveling  in  the  target  room.  The  sabot  can  then  be  allowed  to  impact 
the  target  and  be  destroyed,  usually  without  causing  damage  to  the  target,  or  steel  plates 
can  be  used  to  pulverize  and  deflect  die  sabot  before  it  impacts  the  target 

The  velocity  measurement  system  of  the  range  consists  of  a  Hall  Intervalometer 
System  which  is  located  alongside  the  drift  tube.  It  is  used  to  determine  both  the 
projectile  velocity  and  orientation.  A  pair  of  shadowgrams  from  orthogonal  viewing 
angles  are  recorded  by  a  streak  camera  as  the  projectile  passes  each  of  two  stations  along 
the  drift  tube.  Figure  6.2  is  an  example  of  the  type  of  photograph  which  is  recorded  with 
the  streak  camera.  The  photograph  shows  top  and  side  views  of  the  projectile  at  two 
stations  that  are  0.75  metres  apart.  A  time  reference  is  supplied  by  timing  marks  placed 
along  the  sides  of  the  film  using  a  time  code  generator.  Measurements  at  each  station  are 
used  to  determine  the  velocity  and  yaw  rate  of  the  projectile  in  free  flight  with  accuracies 
of  0.2  percent  and  ±1  degree,  respectively. 

Targets  are  housed  in  a  5.6-m-long  by  4.6-m-wide  by  3.0-m-high  reinforced- 
concrete  enclosed  target  room.  Viewing  ports  on  the  side  and  in  the  ceiling  of  the  room 
permit  high-speed  motion  pictures  of  projectile/target  interaction;  flash  x-ray  heads  and 
breakscreens  can  be  located  freely  within  the  room  for  additional  experiment  diagnostics. 
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A  penetration  experiment  typically  involves  the  launch  of  a  steel  projectile  into  a 
target  consisting  of  geologic  and/or  man-made  materials  contained  within  steel  culverts. 
Projectiles  are  fabricated  in  the  WES  Machine  Shop  by  machining  the  projectile  to  a 
slight  oversize  dimension  from  steel  stock,  heat  treating  it  to  the  desired  hardness  and 
then  remachining  it  to  final  dimensions.  Complexity  of  the  target  can  range  from  a 
simulated  half-space  of  concrete,  to  a  layered  system  of  concrete  and  soil,  to  a  subscale 
model  of  an  actual  structure  or  components  of  the  structure.  Targets  can  be  fabricated  by 
placing  concrete  that  has  been  mixed  by  WES  personnel,  mixed  at  a  local  batch  plant 
under  supervision  of  WES  personnel,  or  simply  purchased  from  a  local  vendor;  soil  fills 
can  be  placed  to  desired  specifications. 

Basic  information  obtained  during  projectile  penetration  experiments  includes 
impact  velocity  and  projectile  orientation  from  the  Hall  Station,  depth  of  penetration,  and 
target  damage  such  as  crater  profile,  crack  patterns  and  photographs.  Additional 
information  can  include  high  speed  movies  (Figure  6.3.a),  flash  x-rays  (Figure  6.3.b),  and 
instrumentation  of  the  projectile  using  an  accelerometer  and  miniature  hardened  data 
acquisition  package  to  record  the  loading  history  (Figure  6.3. c)  which  can  be  integrated  to 
determine  the  projectile’s  motion-time  response  during  the  penetration  process.  Concrete 
cylinders  poured  during  target  fabrication  are  broken  on  or  within  a  few  days  of  the 
impact  test  to  determine  the  rmconfined  compressive  strength  of  the  concrete.  Additional 
experiments  to  determined  the  mechanical  properties  of  the  concrete  can  be  conducted  on 
specimens  cored  from  samples  of  the  concrete  that  were  obtained  during  target 
fabrication. 

6.1.2  Penetration  Experiments 

A  series  of  penetration  experiments  were  conducted  by  launching  AP  projectiles 
into  targets  consisting  of  CSPC  concrete  at  impact  velocities  (V;)  ranging  from  277  to  800 
m/s.  All  experiments  were  conducted  at  normal  incidence  to  the  target  face.  The  CSPC 
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concrete  targets  were  cast  in  corrugated  steel  culverts  approximately  1.37  m  or  1.22  m  in 
diameter.  The  diameter  of  the  target  should  be  greater  than  25  projectile  diameters 
(calibers)  to  eliminate  edge  effects  during  the  penetration  process.  Target  diameters  for 
these  experiments  were  a  minimum  of  45  calibers.  The  length  of  the  targets  varied 
depending  on  impact  velocity,  but  all  target  lengths  were  greater  than  twice  the  depth  of 
penetration  to  prevent  backface  effects  on  the  projectile  penetration.  Cylinders  measuring 
152  mm  in  diameter  by  305  mm  in  length  were  cast  from  the  concrete  placed  in  each 
target  for  subsequent  vmconfined  compressive  strength  tests.  Values  of  unconfined 
compressive  strength  for  each  target  are  given  in  Table  6. 1 .  Penetration  experiments  and 
unconfined  compressive  strength  tests  were  conducted  at  about  28  days  after  the  targets 
were  cast.  Target  strength  was  nominally  36  MPa  and  ranged  from  32.4  to  40.1  MPa. 

The  AP  projectiles  were  dimensioned  as  shown  in  Figure  6.4,  and  were  machined 
from  4340  steel  rods  and  heat  treated  to  a  Rockwell  hardness,  R,.,  of  43-45.  The  iimer 
cavity  of  the  projectile  was  filled  with  either  grout  or  sand  at  a  density  of  approximately 
1.58  Mg/m^  The  projectiles  were  fitted  with  plastic  sabots  and  obturators  that  fit  snugly 
into  the  gun  bore.  The  sabots  and  obturators  were  stripped  aerodynamically  prior  to 
impacting  the  targets.  Pitch  and  yaw  angles  measured  from  the  Hall  Intervelometer  streak 
film  were  less  than  one  degree.  Data  from  the  experiments  is  presented  in  Table  6.1. 

The  primary  data  from  these  experiments  are  depth  of  penetration  and  a  mapping  of 
the  impact  crater  and  penetration  path.  Except  for  the  experiment  at  277  m/s,  all 
projectiles  penetrated  at  least  one  projectile  body  length  into  the  target.  Experiment 
results  are  presented  as  depth  of  penetration  in  calibers  (depth/projectile  diameter)  versus 
impact  velocity  in  Figure  6.5.  The  data  show  a  consistent  trend,  as  expected,  of  increased 
depth  of  penetration  with  increasing  impact  velocity.  The  mapping  of  the  impact  crater 
and  penetration  path  for  the  experiments  conducted  at  impact  velocities  of 277, 499,  and 
642  m/s  is  shown  in  Figure  6.6.  The  impact  crater  enters  the  tunneling  phase  after  a 
penetration  of  about  50  mm,  or  about  2  projectile  diameters.  Posttest  photographs  of  the 
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Table  6.1.  Summary  of  penetration  experiments  for  the  AP  projectile  into  CSPC  concrete. 
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Figure  6.5.  Depth  of  penetration  for  AP  projectiles  impacting  CSPC  concrete. 
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targets  (Figure  6.7)  illustrate  that  the  craters  are  similar  in  appearance  for  all  impact 
velocities,  but  the  number  of  visible  radial  cracks  increases  as  impact  velocity  increases. 
The  cracks  occxir  after  the  penetration  event  and  do  not  influence  the  final  depth  of 
penetration. 

The  targets  for  the  six  experiments  above  were  opened  and  the  projectiles  were 
removed  to  inspect  their  condition.  Posttest  photographs  of  the  six  projectiles  are  shown 
in  Figure  6.8.  Target  material  on  the  front  half  of  the  projectiles  is  believed  to  have 
adhered  to  the  projectile  after  it  had  stopped.  All  of  the  projectiles  show  very  little 
damage  so  that  the  projectiles  can  be  modeled  as  being  effectively  rigid. 

6.1.3  Perforation  Experiments 

The  targets  used  in  the  perforation  experiments  consisted  of  WES5000  concrete 
having  a  nominal  unconfined  compressive  strength  of  38.2  MPa.  The  unreinforced 
concrete  slabs  had  thicknesses  of 284, 254, 216  and  127  mm.  The  concrete  was  poured 
into  steel  culverts  having  a  nominal  diameter  of  1 .52  m,  vibrated  to  remove  trapped  air, 
and  allowed  to  cure  approximately  28  days  prior  to  the  perforation  experiment.  Water 
was  ponded  on  the  top  surface  of  the  slabs  for  approximately  7  days.  The  target  diameter 
is  30  times  the  projectile  diameter  so  that  edge  effects  should  not  influence  the 
experiment  results. 

A  drawing  and  the  pertinent  characteristics  of  the  projectile  used  in  the  perforation 
experiments  is  presented  in  Figure  6.9.  The  projectile  was  machined  fi*om  4340  steel  rods 
and  then  heat  treated  to  a  Rockwell  hardness,  R^jOf 43-45.  The  itmer  cavity  of  the 
projectiles  was  filled  with  sand.  The  masses  of  the  test  projectiles  were  similar  (see  Table 
6.2).  Impact  velocity  (Vj)  for  all  experiments  was  about  313  m/s. 


591  m/s  e.  Vi  =  642m/s  f.  Vj  =  800m/s 

Figure  6.7.  Posttest  photographs  of  CSPC  concrete  targets  penetrated  by  AP  projectiles. 
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Primary  data  obtained  from  these  experiments  included  high-speed  movies  of  the 
slab  impact  viewed  from  the  side  at  an  angle  of  about  45  degrees,  high-speed  movies  of 
the  damage  to  the  backface  of  the  slab  viewed  using  a  mirror  placed  at  an  angle  of  45 
degrees  to  the  backface,  exit  velocity  (Vg)  from  breakscreens  and  the  high-speed  movies, 
and  mappings  of  the  impact  and  exit  craters.  All  experiments  were  conducted  at  normal 
incidence  to  the  slab.  From  the  high-speed  movies  of  the  impact  face,  the  yaw  at  impact 
in  the  vertical  plane  was  less  than  1 .3  degrees  for  all  experiments  except  PERF-7  (2.1 
degrees);  yaw  in  the  horizontal  plane  at  impact  is  not  known  but  is  believed  to  be  small. 

A  typical  target  configuration  is  shown  in  Figure  6.10.  Figme  6.10.a  shows  the  frontface 
of  the  target  with  a  break  screen  attached  to  determine  the  time  of  impact.  Figure  6.  lO.b 
shows  the  backface  of  the  target  with  a  152-mm-square  grid  painted  on  the  surface  to  aid 
in  viewing  the  breakup  of  the  target  as  it  is  being  perforated.  Figures  6. 1 0.c  and  6. 1 0.d 
show  the  posttest  condition  of  the  target.  The  impact  crater  is  significantly  smaller  than 
the  exit  crater.  Also,  most  of  the  cracks  that  were  foimd  on  the  front  of  the  target  were 
also  found  on  the  back  of  the  target  indicating  that  they  extended  through  the  target. 

Projectile  masses  and  impact  and  exit  velocities  are  summarized  in  Table  6.2.  In 
all  of  the  experiments  the  projectile  exited  the  backface  of  the  slab.  For  experiments  on 
the  284-mm-thick  slab,  the  projectile  was  found  on  the  floor  just  behind  the  slab  and  is 
believed  to  have  simply  fallen  through  the  back  of  the  slab  due  to  gravity. 

A  plot  of  VgA^i  versus  slab  thickness  (T)  is  shown  in  Figure  6.1 1 .  Assuming  a 
constant  deceleration  while  perforating  the  slabs,  the  thickness  at  which  the  projectile  will 
just  exit  the  target  is  about  260  mm.  Selected  images  captured  from  the  high-speed 
movies  of  the  impact  face  of  a  target  are  shown  in  Figure  6.12.  From  these  images,  only  a 
cloud  of  "pulverized"  material  and  small  pieces  of  concrete  can  be  seen  coming  from  the 
impact  area.  Selected  images  from  the  high-speed  movies  of  the  backface  of  the  slab  are 
shown  in  Figures  6.13  and  6.14  for  254-  and  127-mm-thick  slabs,  respectively.  These 
images  indicate  that  damage  begins  with  cracks  radiating  from  the  center  of  the  backface. 
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.  Posttest  photo  of  frontface  d.  Posttest  photo  of  backface 

Figure  6.10.  Pre-  and  posttest  phototgraphs  from  a  typical  perforation  experiment. 
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Figure  6.11.  Ratio  of  exit  velocity  to  impact  velocity  versus  target  thickness  for  a  SAP 
projectile  perforating  WES5000  concrete  slabs. 
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C.  t3>t2  d.  t4>t3 

Figure  6.12.  Selected  images  showing  the  front  of  a  254-nim-thick  WES5000  concrete  slab  being  impacted  by  a  SAP  projectile. 


C.  tj  >  t2  d.  t4  >  tj 

Figure  6.13.  Selected  images  showing  the  back  of  a  254-mm-thick  WES5000  concrete  slab  being  perforated  by  a  SAP  projectile. 


.  Selected  images  showing  the  back  of  a  127-mm-thick  WES5000  concrete  slab  being  perforated  by 
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The  outer  edges  of  the  exit  crater  begin  to  form  as  material  is  pushed  out  from  the 
backface  of  the  slab.  As  the  tip  of  the  projectile  protrudes  through  the  backface  of  the 
slab,  the  material  ahead  of  the  projectile  is  destroyed.  The  pieces  being  pushed  from  the 
backface  of  the  slab  were  generally  larger  and  fewer  in  number  for  the  thicker  slabs. 
Pieces  of  concrete  seen  in  the  movies  and  recovered  after  the  experiment  were  relatively 
thick  and  showed  no  delamination  of  concrete  near  the  backface  that  might  result  from  a 
reflected  tensile  stress  wave.  Profiles  showing  the  impact  and  exit  craters  for  three  slab 
thicknesses  are  shown  in  Figure  6. 1 5.  The  size  of  the  exit  crater  increases  with  increasing 
thickness  of  the  slab.  Several  of  the  projectiles  exited  the  slabs  with  significant  yaw  rates. 
Since  the  yaw  at  impact  was  small,  the  exit  yaw  rates  were  attributed  to  inhomogeneity  of 
the  concrete  and  interaction  of  the  projectile  with  the  pieces  of  concrete  being  pushed 
through  the  backface  of  the  slab. 

6.2  NUMERICAL  SIMULATIONS 

Simulations  of  the  penetration  and  perforation  experiments  were  conducted  using 
the  NIF  model  implemented  in  the  EPIC  code.  Initial  simulations  of  the  penetration 
experiments  indicated  some  instability  that  was  traced  to  the  hydrostatic  tension  and 
unloading/reloading  portions  of  the  model.  In  Figure  4.1,  the  hydrostatic  tension  response 
is  nonlinear  as  the  pressure  is  approached.  This  response  was  altered  to  be  linear  with 
a  slope  equal  to  as  defined  in  Section  4.2. 1  and  a  pressure  limit  equal  to  The 
shear  unload/reload  portion  was  altered  to  be  linear  with  a  slope  equal  to  1 .20  times  the 
initial  loading  slope  of  the  corresponding  Yoc/  curve.  All  simulations  were 
axisyimnetric  and  used  triangle  elements.  The  simulations  were  performed  on  the  WES 
Cray  C90  super  computer.  Post-processing  was  performed  using  software  supplied  with 
EPIC  and  downloaded  as  encapsulated  post-script  files  for  insertion  into  the  figures. 
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a.  127-mm-thick  slab  b.  216-mm-thick  slab  c.  254-mm-thick  slab 

Figure  6.15.  Impact  (left)  and  exit  (right)  crater  profiles  for  the  perforation  experiments. 
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6.2.1  Penetration  Experiments 

The  AP  projectile  penetration  experiments  into  CSPC  concrete  conducted  at  impact 
velocities  of 211, 499,  and  642  m/s  will  be  simulated  using  the  NIF  model.  A  close-up  of 
the  initial  grid  used  in  the  simulations  is  shown  in  Figure  6.16.a.  Figure  6.16.b  shows  the 
materials  used  in  the  grid.  In  the  projectile,  the  outer  two  columns  of  the  grid  make  up 
the  steel  case  of  the  projectile.  Although  not  shown,  a  steel  end-cap  sealed  the  back  of 
the  projectile.  The  sliding  interface  algorithm  used  in  the  code  to  separate  the  projectile 
case  from  the  concrete  target  requires  that  the  target  elements  be  about  the  same  size  or 
smaller  than  the  outer  elements  of  the  projectile.  This  aids  in  ensuring  that  target 
elements  do  not  “mix”  with  the  projectile  elements.  Since  the  projectile  case  is  relatively 
thick  and  no  significant  deformation  is  expected,  the  influence  of  movement  of  the  fill 
material  on  the  projectile  response  should  be  minor.  The  fill  material  was  “fixed”  to  the 
steel  case  so  that  no  slideline  was  required. 

The  full  initial  grids  for  the  simulations  are  shown  in  Figure  6.17.  The  same  grid 
was  used  for  both  the  277  and  499  m/s  simulations.  Beyond  a  radius  of  102  mm  (4 
inches)  the  grid  was  gradually  expanded  to  a  final  radius  of  508  mm  (20  inches). 

Similarly,  beyond  a  depth  of  559  mm  (22  inches)  the  grid  was  gradually  expanded  to  a 
final  depth  of  1,016  mm  (40  inches).  The  lateral  extent  of  the  grid  for  the  642  m/s 
simulation  was  the  same  as  that  for  the  two  slower  impact  velocity  simulations.  Since  the 
depth  of  penetration  is  expected  to  be  greater  for  the  higher  impact  velocity,  the  total 
depth  of  the  grid  was  extended  to  2,032  mm  (80  inches)  with  the  expanding  grid 
beginning  at  a  depth  of  1 , 1 1 8  mm  (44  inches).  The  grid  was  expanded  to  reduce  the  size 
of  the  model  and  to  decrease  computation  time.  Element  sizes  were  kept  relatively  small 
within  the  region  of  the  target  influencing  the  penetration  of  the  projectile. 

Model  parameters  presented  in  Table  4.2  were  used  in  the  simulation.  The  value  of 
DYN  applied  to  the  ultimate  strength  in  both  shear  and  hydrostatic  tension  was  1 .25.  A 
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model  parameter  that  is  not  defined  by  mechanical  property  tests  is  the  critical  value  of 
the  strain  ratio  EBAR  that  was  discussed  in  Section  5. 1 .2.  This  value  must  initially  be 
determined  for  a  given  target  material  by  comparing  a  simulation  result  with  an 
experiment  result.  Currently,  this  must  be  done  through  trial-and-error.  The  critical 
value  of  EBAR  (EBAR^rit)  was  determined  to  be  400  for  the  CSPC  concrete  by  simulating 
the  penetration  experiment  at  277  m/s  imtil  a  reasonable  result  (depth  of  penetration)  was 
obtained.  This  value  of  EBARj^t  was  then  used  for  the  other  two  impact  velocity 
simulations.  During  initial  simulations  with  the  model,  small  “sliver”  elements  that  had 
detached  from  the  grid  were  overlapping  with  the  remaining  target  mesh.  Although  these 
elements  were  not  influencing  the  simulation  results,  they  were  causing  the  time  step  to 
be  uimecessarily  small.  The  strain  in  these  elements  was  not  sufficient  to  allow  normal 
removal  using  EBAR^rit  criterion.  In  order  to  gradually  remove  these  elements,  each 
element  was  checked  to  see  if  its  minimum  height  was  less  than  1.27  mm  (0.05  inches). 

If  this  check  was  true,  the  value  of  EBAR  for  that  element  was  multiplied  by  100.  If  the 
increased  value  of  EBAR  was  greater  than  EBARcnt  then  the  element  was  removed  from 
the  simulation.  This  increased  computation  performance  without  significantly  affecting 
the  simulation  results. 

Depth  of  penetration  from  the  simulations  is  compared  to  the  experiment  results  in 
Figure  6.18.  Both  the  simulations  and  the  experiments  show  an  increase  in  depth  of 
penetration  as  impact  velocity  increases.  The  simulations  resulted  in  slightly  greater 
depth  of  penetration  than  the  experiments.  The  depths  from  the  simulations  are 
approximately  20  percent  greater  than  the  depths  from  the  experiments  at  impact 
velocities  of  277  and  499  m/s  and  40  percent  greater  at  642  m/s.  It  is  possible  that  the 
simulation  results  could  be  improved  by  modifying  EB AR^nt  to  better  fit  the  experiment 
results  at  one  of  the  other  impact  velocities. 

The  deformed  grid  at  several  times  during  the  simulation  is  shown  in  Figures  6.19, 
6.20,  and  6.21  for  the  277, 499,  and  642  m/s  impact  velocities,  respectively.  The  images 
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Figure  6.18.  Comparison  of  simulation  and  experiment  results  for  an  AP  projectile  impacting  CSPC  concrete. 
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Figure  6.19.  Deformed  grid  for  an  AP  projectile  impacting  CSPC  concrete  at  277  m/s 
(time  is  in  seconds;  axis  units  are  inches). 
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Figure  6.20.  Deformed  grid  for  an  AP  projectile  impacting  CSPC  concrete  at  499  m/s 
(time  is  in  seconds;  axis  units  are  inches). 
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Figure  6.21.  Deformed  grid  for  an  AP  projectile  impacting  CSPC  concrete  at  642  m/s 
(time  is  in  seconds;  axis  units  are  inches). 
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show  the  gradual  formation  of  the  impact  crater  and  penetration  tunnel.  As  the 
penetration  progresses,  material  around  the  impact  area  is  ejected  in  a  manner  similar  to 
that  shown  in  Figures  2.3  and  6.12.  The  diameter  of  the  impact  crater  at  the  target  surface 
for  the  277  m/s  impact  velocity  is  about  200  mm  (7.9  inches).  The  crater  being  generated 
in  the  simulation  shown  in  Figure  6. 19  is  about  254  mm  (10  inches)  in  diameter.  After 
the  impact  phase,  both  the  simulation  and  experiment  show  a  tunnel  that  is  the  same 
diameter  as  the  projectile.  The  deformed  grids  in  Figures  6.20  and  6.21,  however,  exhibit 
continued  deformation  along  the  penetration  path  after  the  projectile  has  passed.  This 
results  in  a  larger  impact  crater  and  tunnel  than  was  seen  in  the  experiments.  The  cause 
of  this  difference  is  the  criteria  used  in  the  simulations  to  delete  highly  damaged  and  very 
small  elements  from  the  finite  element  grid  in  order  to  improve  computational  time.  In 
the  simulations  the  extent  of  damage  and  degradation  in  the  target  is  indicated  by  the 
value  of  the  parameter  EBAR  (strain  ratio  in  the  NIF  model).  Figures  6.22, 6.23,  and 
6.24  show  the  progression  of  the  change  in  EBAR  at  two  times  during  the  simulations  for 
each  of  the  impact  velocities.  The  simulations  show  high  values  of  EBAR  indicating 
significant  damage  and  degradation  of  the  material  where  the  impact  crater  and  tunnel  are 
being  formed  and  aroimd  the  aft  end  of  the  projectile.  The  radial  extent  of  the  material 
with  high  values  of  EBAR  is  about  four  times  the  diameter  of  the  projectile.  Beyond  this 
range,  the  material  is  effectively  intact  and  fairly  competent.  In  the  simulations,  the 
material  that  is  highly  damaged  or  degraded  is  eventually  removed  from  the  grid  if  EBAR 
for  the  element  exceeds  EB  ARgn,.  The  impact  crater  and  tunnel  continue  to  grow  with 
removal  of  the  damaged  elements  and  eventually  become  larger  than  the  impact  crater 
and  tunnel  observed  in  the  experiments.  In  the  experiments,  material  along  the  tunnel  is 
highly  damaged  but  not  ejected  from  the  target  because  it  is  held  in  place  by  its 
neighboring  material  due  to  interlocking  and  fiictional  forces.  During  removal  of 
projectiles  from  the  targets  after  the  penetration  experiments,  it  was  noted  that  material 
along  the  penetration  path  was  very  weak  and  easily  broken. 
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6.2.2  Perforation  Experiments 

Perforation  experiments  conducted  into  the  254-,  21 6-,  and  127-nim-thick  slabs 
will  be  simulated  using  the  NIF  for  the  target  response.  A  close-up  of  the  initial  grid  used 
in  the  simulations  of  the  perforation  experiments  is  shown  in  Figure  6.25. a.  Figure  6.25.b 
shows  the  materials  used  in  the  grid.  In  the  projectile,  the  outer  column  of  the  grid  makes 
up  the  steel  case  of  the  projectile.  Although  not  shown,  a  steel  end-cap  sealed  the  back  of 
the  projectile.  No  significant  deformation  of  the  projectile  case  is  expected,  so  the 
influence  of  movement  of  the  fill  material  on  the  projectile  response  should  be  minor. 

The  fill  material  was  “fixed”  to  the  steel  case  so  that  no  slideline  was  required. 

Model  parameters  presented  in  Table  4.4  were  used  in  the  simulations.  The  value 
of  DYN  applied  to  the  ultimate  strength  in  both  shear  and  hydrostatic  tension  was  1 .25. 

In  order  to  determine  the  value  of  EBARcrit,  a  simulation  of  the  SAP  projectile  impacting 
a  half-space  of  WES5000  concrete  was  performed.  Experience  has  shown  that  the 
projectile  should  penetrate  approximately  3.7  calibers  into  this  type  of  concrete.  The 
value  of  EBAR^ri,  was  determined  to  be  350  by  simulating  a  penetration  experiment  at 
313  m/s  xmtil  a  reasonable  result  (depth  of  penetration)  was  obtained.  The  deformed  grid 
firom  the  final  simulation  is  shown  in  Figure  6.26.  The  final  depth  of  penetration  is  about 
3.77  calibers.  This  value  was  then  used  for  the  three  simulations  of  the  perforation 
experiments.  Modification  to  the  tension  and  shear  unload/reload  responses  and  the 
additional  check  to  gradually  remove  elements  with  minimum  height  less  than  1.27  mm 
(0.05  inches)  discussed  above  were  also  used  in  these  simulations. 

The  initial  grids  for  the  perforation  simulations  are  shown  in  Figure  6.27.  The 
radius  of  the  grids  were  all  the  same  and  equal  to  762  mm  (30  inches).  Since  the  size  of 
the  targets  were  relatively  small  it  was  not  necessary  to  expand  the  grids.  All  elements 
within  the  grids  were  the  same  size.  The  impact  velocity  for  all  simulations  was  313  m/s. 
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Figure  6.25.  Close-up  of  initial  grid  and  materials  used  in  simulation  of  the  SAP  projectile  perforating  WES5000  concrete  slabs 
(axis  units  are  inches). 
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Figure  6.26.  Deformed  grid  for  a  SAP  projectile  impacting  a  half-space  of  WES5000 
concrete  (time  is  in  seconds;  axis  units  are  inches). 
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Figure  6.27.  Initial  grids  for  simulating  the  SAP  projectile  perforating  127-,  216-,  and 
254-nim-thick  (5.0-,  8.5-,  and  10.0-inch-thick)  WES5000  concrete  slabs 
(axis  units  are  inches). 
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.  Each  of  the  simulations  exhibited  a  trend  of  decreasing  velocity  with  time  until  a 
constant  exit  velocity  was  achieved.  Results  from  the  simulations  are  compared  to  the 
experiment  results  in  Figure  6.28  by  comparing  the  ratio  of  exit  velocity  to  the  impact 
velocity.  The  simulation  results  agree  well  with  the  experiment  results,  with  the  greatest 
difference  being  for  the  254-mm-thick  slab.  The  deformed  grids  at  several  times  during 
the  simulations  for  each  of  the  slab  thicknesses  are  shown  in  Figures  6.29, 6.30,  and  6.31, 
respectively.  Contour  plots  of  EB  AR  (strain  ratio)  for  each  of  the  target  thicknesses  are 
shown  in  Figures  6.32,  6.33,  and  6.34,  respectively.  Like  the  experiment  results  in  Figure 
6.12,  each  simulation  shows  the  formation  of  the  impact  crater  as  material  is  ejected  from 
the  impact  area.  As  the  thickness  of  the  slab  increases,  the  size  of  the  exit  crater  increases 
as  shown  in  Figures  6.13, 6.14,  and  6.15.  Simulations  for  the  216-  and  254-mm-thick 
slabs  also  show  the  progressive  formation  of  the  exit  cone  as  was  described  in  Section 
2.4.  As  the  simulation  progresses,  the  projectile  passes  through  the  cone  essentially 
destroying  it.  Large  areas  of  “damaged”  material  around  the  impact  crater,  projectile,  and 
exit  crater  are  implied  by  the  strain  ratio  contours  in  Figures  6.32, 6.33,  and  6.34.  As 
with  the  penetration  experiments,  material  within  these  regions  was  easily  broken,  but 
farther  away,  the  material  was  very  competent.  Again,  the  highly  “damaged”  material  in 
the  simulations  continues  to  be  removed  if  the  value  of  EBAR  exceeds  EBARerit-  The 
continued  removal  of  these  highly  “damaged”  elements  allows  the  impact  and  exit  craters 
to  continue  to  grow. 
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Figure  6.29.  Deformed  grid  for  a  SAP  projectile  perforating  a  127-mm-thick  WES5000 
concrete  slab  (time  is  in  seconds;  axis  units  are  inches). 
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Figure  6.33.  EBAR  (strain  ratio)  contours  for  a  SAP  projectile  perforating  a  216-mm- 
thick  WES5000  concrete  slab  (time  is  in  seconds;  axis  units  are  inches; 
contours  are  dimensionless). 
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Figure  6.34.  EBAR  (strain  ratio)  contours  for  a  SAP  projectile  perforating  a  254-mm- 
thick  WES5000  concrete  slab  (time  is  in  seconds;  axis  units  are  inches; 
contours  are  dimensionless). 
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CHAPTER  7 

CONCLUSIONS  AND  RECOMMENDATIONS 

7.1  CONCLUSIONS 

Simulation  of  high-velocity  projectile  penetration  into  brittle  geomaterials  is 
generally  performed  using  either  empirical,  analytical,  or  numerical  techniques.  Only  the 
numerical  techniques  offer  the  capability  of  detailed  simulation  of  the  phenomena  that 
occur  including  the  large  deformation  of  the  target  during  formation  of  craters  and  the 
tunneling  phase,  as  well  as  insight  into  the  response  and  condition  of  the  target  as  a  result 
of  the  penetration  event.  The  loads  applied  to  the  target  during  these  events  are  intense 
with  pressures  approaching  1,000  MPa  and  strain  rates  up  tolO'*  /sec. 

Laboratory  mechanical  property  tests  on  brittle  geomaterials  show  a  transition  in 
material  response  from  brittle  at  low  pressures  to  ductile  at  high  pressures.  A  nonlinear, 
inelastic  fracture  (NIF)  model  was  developed  that  can  replicate  the  brittle  response  of  the 
material  at  low  pressures  and  the  transition  to  ductile  response  at  high  pressures.  The 
underlying  hypothesis  of  the  model  is  that  the  shearing  response  of  the  material  can  be 
resolved  into  a  brittle  cohesive  component  and  a  ductile  frictional  component.  At  very 
low  pressures  the  cohesive  component  controls  the  behavior  of  the  material  and  the 
response  is  brittle.  As  pressure  increases,  the  cementing  bonds  between  the  aggregate 
particles  are  broken  and  the  contribution  of  the  cohesive  component  decreases  while  the 
contribution  of  the  frictional  component  increases.  This  phase  is  the  transition  from 
brittle  to  ductile  response.  The  response  becomes  fully  frictional  once  the  cementing 
bonds  are  completely  broken.  This  is  the  fully  ductile  phase  of  the  response.  The  NIF 
model  response  agreed  well  with  the  results  from  various  quasi-static  triaxial  experiments 
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on  concrete  samples  and  captured  all  three  of  the  shear  phases.  Dynamic  loading 
experiments  indicate  strength  enhancement  by  factors  of  two  or  more  as  loading  rate 
increases.  Dynamic  strain  rate  effects  were  indirectly  accoimted  for  in  the  model  by 
multiplying  the  ultimate  strength  of  the  material  in  shear  and  hydrostatic  tension  by  a 
constant  parameter. 

The  model  was  implemented  into  a  finite-element  wave  propagation  code.  This 
process  involved  implementing  algorithms  to  calculate  the  correct  strain  input  to  the 
model,  and  to  convert  the  stress  output  by  the  model  into  the  form  accepted  by  the  code. 
The  model  requires  the  total  Lagrangian  strain  based  on  the  original  configuration  as 
input.  The  Green-Lagrange  measure  of  strain  meets  this  requirement  and  was  selected  as 
the  strain  input.  The  strain  components  were  calculated  using  the  deformation-rate  and 
spin  tensors  fi-om  the  updated  Lagrangian  kinematics  formulation  of  the  code.  The  stress 
output  by  the  model  is  the  Cauchy  stress  aligned  with  the  original  configuration.  This 
stress  was  rotated  to  the  current  configuration  of  the  element  to  be  compatible  with  the 
formulation  of  the  code. 

A  series  of  laboratory  penetration  and  perforation  experiments  were  conducted  to 
evaluate  the  performance  of  the  model  within  the  finite-element  code.  Penetration 
experiments  were  conducted  by  laimching  a  robust,  thick-walled  steel  projectile  into 
semi-infinite  concrete  targets  at  impact  velocities  ranging  from  277  to  800  m/s.  The  data 
from  these  experiments  included  depth  of  penetration  and  a  mapping  of  the  impact  crater 
and  penetration  path.  Perforation  experiments  were  conducted  by  launching  a  robust  steel 
projectile  at  3 13  m/s  into  concrete  slabs  with  thicknesses  ranging  from  127  to  284  mm. 
Data  from  these  experiments  included  high-speed  movies  of  the  impact  phase,  high-speed 
movies  to  capture  the  evolution  of  damage  to  the  backface  of  the  slabs  and  to  measme  the 
exit  velocity,  and  the  mappings  of  the  impact  and  exit  craters. 
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Simulation  of  penetration  experiments  into  thick  concrete  targets  showed  the 
formation  of  the  impact  crater  and  the  tuimel  during  the  penetration  event.  The  depth  of 
penetration  from  the  simulations  agreed  well  with  the  data  from  the  experiments.  The 
simulations  showed  significant  damage  and  degradation  of  the  material  where  the  impact 
crater  and  tunnel  were  being  formed  and  around  the  aft  end  of  the  projectile.  As  the 
simulation  continued,  the  impact  crater  and  tunnel  continued  to  grow  and  eventually 
became  larger  than  the  impact  crater  and  tuimel  observed  from  the  experiments.  During 
posttest  examination  of  the  targets  from  the  experiments,  it  was  noted  that  material  along 
the  penetration  path  was  veiy  weak  and  easily  broken.  The  broken  materials  were  not 
ejected  from  the  target  because  they  were  held  in  place  by  interlocking  and  frictional 
forces,  which  were  not  accounted  for  in  the  simulations.  Exit  velocities  from  simulation 
of  the  perforation  experiments  agreed  well  with  the  data  from  the  experiments.  The 
smudations  showed  the  formation  of  the  impact  crater  and  the  progressive  formation  of 
the  exit  cone  ahead  of  the  projectile  as  the  backface  of  the  target  was  approached. 

7.2  RECOMMENDATIONS 

The  criteria  used  to  eliminate  highly  deformed  and  damaged  elements  from  the 
simulations  need  further  development.  These  criteria  are  needed  to  allow  selection  of 
reasonable  computational  time  increments.  Currently,  determination  of  the  parameters 
controlling  these  criteria  is  iterative  and  requires  comparison  of  simulation  results  with 
data  from  penetration  experiments.  Since  values  for  EB  ARj.rit  were  similar  for  both  of  the 
concrete  materials  used  in  this  research,  a  value  in  the  range  of  350  to  400  may  be 
sufficient,  but  this  requires  further  investigation.  A  criterion  that  does  not  require  this 
iterative  approach  or  can  be  based  purely  on  the  response  and  condition  of  the  damaged 
material  is  preferable.  Use  of  the  model  with  meshless  methods  and  adaptive  grids  may 
eliminate  the  need  for  removing  highly  deformed  elements  and  is  a  topic  for  future 
investigation.  Since  these  methods  can  require  significant  computational  time,  parallel 
processing  may  be  required. 
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Additional  data  on  the  response  of  brittle  geomaterials  to  the  combined  effects  of 
pressure  and  strain  rate  are  needed.  Also,  current  data  is  inconclusive  as  to  whether  strain 
rate  affects  the  failure  mode  of  the  material  or  just  its  ultimate  strength. 

The  ultimate  strength  in  the  current  model  does  not  distinguish  between 
compression  and  extension  loading.  Data  from  triaxial  experiments  on  brittle 
geomaterials  show  that  the  strength  in  extension  is  much  lower  than  the  strength  in 
compression.  Use  of  different  compression  and  extension  strength  surfaces  based  on  the 
Lode  angle  would  enhance  the  capabilities  of  the  model. 

Other  enhancements  to  the  model  can  include  thermal  effects  and  shear  induced 
dilatancy.  Significant  heat  is  generated  diiring  high-velocity  projectile  penetration,  but  it 
is  unclear  as  to  the  effect  heat  may  have  on  the  material  response.  Shear  induced 
dilatancy  can  be  added  to  the  model  by  including  a  correction  function  to  account  for  the 
volumetric  strains  produced  during  deviatoric  loadings. 

Cause  and  effect  studies  can  be  conducted  using  the  model  to  investigate  the 
effect  of  the  cohesive  and  frictional  components  on  the  penetration  process.  Ofiier 
geomaterials  that  are  as  strong  as  concrete  but  not  as  brittle,  such  as  rock,  or  stronger  and 
more  brittle,  such  as  high-strength  concrete,  should  be  fit  to  the  model  and  used  in 
simulations  to  evaluate  the  model  performance  for  these  materials.  Evaluation  of  the 
model  performance  xmder  loading  conditions  resulting  from  explosions  or  quasi-static 
applications  is  also  of  interest. 
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*$ 

SUBROUTINE  NIF (LI , LN, INCOMP, DVBAR, EDEV, EDOT, EXDOT, EXYDOT, 

2  ENSUM, EXZDOT, EYDOT, EYZDOT, EZDOT, SBAR, SMAX, SPINRZ, 

3  SS2 , XSPIN, YSPIN, ZSPIN, EBAR , EP , EPDOT , ES , 

4  ICHECK,M, NODES , NODE4 , SX , S Y , S  Z , SXY , SXZ , S YZ , 

5  BULK,DVOL,U,HMIN,Q,QMAX) 

* 

C  NIF  computes  the  stresses  for  the  JDC  Nonlinear, Inelastic  Fracture 
C  model 

C 

C  called  by  USTRES 

*  Packaged  in  EPIC5 
C 

*  IMPLICIT  NONE 

INTEGER  LI , LN, INCOMP , M, NODES , NODE4 

*  LI  =  first  element  of  block 

*  LN  =  last  element  of  block 

*  M  =  the  material  number 

*  INCOMP  flags  incompressible  materials 
INCLUDE  'VECTMX' 

REAL  DVBAR  {MXLB),EDEV  {MXLB),EDOT  (MXLB), EXDOT  (MXLB) , 

2  EXYDOT (MXLB) , EXZDOT (MXLB) , EYDOT  (MXLB) , EYZDOT (MXLB) , 


3 

EZDOT 

(MXLB), SBAR  (MXLB) , SMAX 

(MXLB) , 

4 

SPINRZ (MXLB) , XSPIN  (MXLB), YSPIN  (MXLB), 

5 

ZSPIN 

(MXLB), ENSUM  (MXLB),U 

(MXLB) ,HMIN 

(MXLB) 

6 

Q 

(MXLB),QMAX  (MXLB) 

EDOT ( I )  =  1 

total  strain  rate 

REAL 

EBAR 

(MXLB) , EP  (MXLB) , EPDOT 

(MXLB) , 

2 

ES 

(MXLB),SX  (MXLB),SY 

(MXLB) ,SZ 

(MXLB) , 

3 

SXY 

(MXLB) , SXZ  (MXLB) , SYZ 

(MXLB) ,BULK 

(MXLB) , 

4 

SS2 

(MXLB) ,DVOL{MXLB) 

*  EBAR (I)  =  plastic  strain 

INTEGER  I CHECK (MXLB) 

INCLUDE  'MATERL' 

INCLUDE  'MATERC 
INCLUDE  'MISC 
INCLUDE  'FILES' 

INCLUDE  'NBSPHC 
C 

REAL  C2D9,DENV,DT2,GDT,GDTI,TDT,TGDT,TGDTI, THIRD, GV,CGV 

*  variables  ending  in  V  are  temporaries  for  vector  quantities 

*  which  are  constant  in  the  current  element  block 
INTEGER  I,J 

*  I  is  an  index  into  vector  arrays 

*  J  is  an  index  into  element  arrays 
REAL  DSX(MXLB) ,DSXY( MXLB) ,DSXZ( MXLB) , 

2  DSY(MXLB) ,DSYZ(MXLB) ,DSZ(MXLB) ,SX1(MXLB) , 

3  SXYl (MXLB) ,SXZ1 (MXLB) , SYl (MXLB) , SZl (MXLB) , SYZl (MXLB) 

PARAMETER  (THIRD  =  1. 0/3.0) 

PARAMETER  (C2D9  =  2. 0/9.0) 

REAL  TW03RD 

PARAMETER  (TW03RD  =  2. 0/3.0) 

Q  ********************************************************************** 
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C  *  Nonlinear,  Inelastic  Fracture  Model  * 

c  *  * 

C  *  DONALD  CARGILE  * 

C  *  USAE  WATERWAYS  EXPERIMENT  STATION  * 

C  *  VICKSBURG,  MISSISSIPPI  * 

C  *  * 

Q  ********************************************************************** 


REAL  VHMDAV(MXLJDC) ,  LDWMDA(3,3),  FDGMDA(3,3),  ROTMDA(3,3), 

1  DRMDA(3,3) ,DRIMDA(3,3) ,DETMDA{3,3) ,ETMDA(3,3) 

REAL  DEMDA ( 6 ) ,  SIGMDA(6) ,  EMDA{6),  EJDC{6),  DEJDC{6) 

REAL  SBMDA{MXLB),  HISJDC(8) 

REAL  EPSKKM , ATYCP , AGAMCP , GMOCTM , FLG JDC , AK , AEPSC , AAH , ABH , ACH , 
AEPSL , AKL , ASIGS , ASIGF , AEPSS , AEPSF , AA, AB , AC , ABO , ABl , AMUC , 
E1E2 , E1E3 , E2E3 , EPS12 , EPS13 , EPS23 , GAMOCT, SIGM, 

EPSKKP, EPSKKO , SIGMSU, DSIGM, RNUM, RDEN, R, A, B, X, GAM, 

FALULT , PCRUSH , PFRIC , ATYC , FACTRF , FAL JDC , FRIC , AGAMC , REDUC , 
HARD, COHSN, DEV, TOCTM, GULRL, RATIO, TAUOCT, GMOCTO , 

AKULRL , GMOD , ANGAMC , DEVEPS , EPSKIW , QIV , Q2V , SMALL , 

EBARP , EBARS , AKLAK , AEPSLC , AAAC , FRICDEN , EPSKKD , 

ATYCUR , ATYCPUR , AGAMCUR , DT2  2 , EFAILV , S IGML , AKEPSC 
REAL  DF(3,3) ,EIN{3,3) ,DEIN(3,3) , LDW2 (3 , 3 ) , LDW3 (3,3) , 

1  DLF (3 , 3 ) , STEMPl (3 , 3 ) , STEMP2 (3,3) 

REAL  EO , AUX , AUXl , AUX2 , SIGMG , DETF , CGVJDC , BLK JDC , EPSKK , EPSKKT , 

1  CMOD JDC , STRFCT 

INTEGER  I IN , JIN , KIN , K , IHYPO 
c 

CDIR$  NOVECTOR 

DENV=  DEN(M) 

TDT  =  2.0*DT 
DT2  =  0.5*DT 
GV  =  G{M) 

CGV  =  4.0*GV/3.0 
GDT  =  GV*DT 
TGDT  =  2.0*GDT 
GDTI  =  1.0 /GDT 
TGDTI  =  1.0 /TGDT 
QIV  =  Ql(M) 

Q2V  =  Q2(M) 

SMALL  =  -EPSLON 
C 
C 

C  ASSIGN  MATERIAL  CONSTANTS 

C 

C  VOLUMETRIC  PART 

C 

STRFCT  =  C18(M) 

AK  =  C0{M) 

AEPSC  =  C1{M) 

AAH  =  C2(M) 

ABH  =  C3(M) 

ACH  =  C4(M) 

AEPSL  =  C5{M) 

AKL  =  C6(M) 

ASIGS  =  C7(M)*STRFCT 
ASIGF  =  C17(M) 

AEPSS  =  C9 (M) *STRFCT 
C  AEPSF  =  C10{M) 

ANGAMC  =4.0 
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AEPSF  =  ANGAMC  *  AEPSS 
AKLAK  =  AKL  -  AK 
AEPSLC  =  AEPSL  -  AEPSC 
AKEPSC  =  AK  *  AEPSC 

SIGML  =  AKEPSC  +  AAH  *  AEPSLC  +  ABH  *  AEPSLC* AEPSLC 
&  +  ACH  *  AEPSLC**3.0 

Calculate  pressure  at  crush  and  pressure  at  start  of 
friction  only 

PCRUSH  =  AKEPSC 

PFRIC  =  AKEPSC  -  AAH  *  ABH/ (3 . 0*ACH) 

1  +  (2.0*ABH*ABH*ABH)/(27.0*ACH*ACH) 

DEVIATORIC  PART 

AA  =  Cll (M) *STRFCT 
AB  =  C12{M) 

AC  =  C13 (M) *STRFCT 
ABO  =  C14(M) 

ABl  =  C15(M) 

AMUC  =  C16(M) 

AAAC  =  AA  -  AC 
EFAILV  =  EFAIL(M) 

IF { IGEOM . LE . 7 ) THEN 
1  =  0 

DIR$  VECTOR 
VD$  VECTOR 
DO  15,  J=L1,LN 
1  =  1  +  1 

SBARd)  =  (SX{I)+SZ(I)+SY(I)  )  *  THIRD 

SXl(I)  =  SX(I)  -  SBAR(I) 

SZl(I)  =  SZ{I)  -  SBAR(I)  ' 

SYl(I)  =  SY(I)  -  SBAR(I) 

SXZl(I)  =  SXZ{I) 

SXYl(I)  =  SXY{I) 

SYZl(I)  =  SYZd) 

DSXd)  =  -SXZld)*SPINRZd)*TDT 

DSY (I )  =0.0 

DSZd)  =  -DSXd) 

DSXYd)  =0.0 

DSXZd)  =  (SXld)-SZld)  )*SPINRZd)*DT 
DSYZ (I )  =  0.0 

EDOTd)  =  AMAXKEDOTd)  ,0.0001) 

15  CONTINUE 

C  enddo 

ELSE 

*  here  IGEOM. EQ. 8 

1  =  0 

DO  16,  J=L1,LN 
1  =  1  +  1 

SBARd)  =  (SXd)+SYd)+SZd)  )*THIRD 
SXld)  =  SXd)  -  SBARd) 

SYld)  =  SYd)  -  SBARd) 

SZld)  =  SZd)  -  SBARd) 

SXYld)  =  SXYd) 

sxzid)  =  sxzd) 

SYZld)  =  SYZd) 


nnn  oo  oono 


190 


DSX(I)  =  {YSPIN(I)*SXZ1(I)  -  ZSPIN(I)*SXY1(I) )*TDT 
DSY(I)  =  (ZSPIN(I)*SXY1{I)  -  XSPIN(I) *SYZ1 (I) ) *TDT 
DSZ(I)  =  {XSPIN{I)*SYZ1(I)  -  YSPIN(I)*SXZ1(I) )*TDT 
DSXY(I)  =  {ZSPIN(I)*(SX1(I)-SY1{I) )  +  YSPIN{I) *SYZ1 (I) 
1  -  XSPIN(I)*SXZ1{I) )*DT 

DSXZ(I)  =  {YSPIN{I)*(SZ1(I)-SX1(I) )  +  XSPIN{I) *SXY1 (I) 
1  -  ZSPIN(I)*SYZ1(I) )*DT 

DSYZ(l)  =  (XSPIN(I)*(SY1(I)-SZ1{I) )  +  ZSPIN(I) *SXZl (I) 
1  -  YSPIN{I)*SXY1(I) )*DT 

EDOT(I)  =  AMAXl{EDOT(I) ,0.0001) 

16  CONTINUE 

C  enddo 

ENDIF 

DIR$  NOVECTOR 


1=0 

DO  100  J=L1,LN 
1=1+1 

SS2(I)  =  0.0 

IFdCHECKd)  .LT.O)  GOTO  100 

DO  110  K=1,MXLJDC,1 
VHMDAV ( K ) = VHMDA ( J , K ) 

CONTINUE 

L=D+W 

IF ( IGEOM . LE . 7 ) THEN 
C 
C 

LDWMDAd,l)  =  (EXDOTd)+ENSUMd)  ) 

LDWMDA (1 , 2 ) = ( EXYDOT (I ) * 0 . 5 ) 

LDWMDAd,3)  =  {EXZDOTd)*0.5)  -  SPINRZd) 
LDWMDA (2 , 1 )= (EXYDOT (1) *0 . 5) 

LDWMDA(2,2)  =  (EYDOTd)+ENSUMd)  ) 
LDWMDA(2,3)  =  (EYZDOTd)*0.5) 

LDWMDA(3,l)  =  {EXZDOTd)*0.5)  +  SPINRZd) 
LDWMDA(3 , 2  )  =  (EYZDOTd)  *0 . 5) 

LDWMDA(3,3)  =  (EZDOTd)+ENSUMd)  ) 
c 

DETMDA ( 1 , 1 ) = ( EXDOT (I ) +ENSUM (I ) ) *DT 
DETMDA (1 , 2 ) = ( EXYDOT (I ) * 0 . 5 ) *DT 
DETMDA (1 , 3 )  =  { EXZDOT (I ) *  0 . 5 ) *DT 
DETMDA(2, 1) = (EXYDOT (1) *0 . 5) *DT 
DETMDA ( 2 , 2 ) = ( EYDOT (I ) +ENSUM ( I ) ) *DT 
DETMDA ( 2 , 3 )  =  ( EYZDOT (I ) *  0 . 5 ) *DT 
DETMDA ( 3 , 1 )=( EXZDOT (I )* 0 . 5 ) *DT 
DETMDA ( 3 , 2 ) = ( EYZDOT (I ) * 0 . 5 ) *DT 
DETMDA(3,3)=(EZDOT(I)+ENSUM(I) ) *DT 
c 

DRMDA(1,1)=  1.0 

DRMDA(1,2)=  0.0 

DRMDA(1,3)=  -0.5*SPINRZ(I)*DT 

DRMDA(2,1)=  0.0 

DRMDA(2,2)=  1.0 

DRMDA(2,3)=  0.0 
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c 


C 

c 


c 


c 


c 


DRMDA(3,1)=  0.5*SPINRZ(I)*DT 
DRMDA(3,2)=  0.0 
DRMDA{3,3)=  1.0 

DRIMDA{1,1)=  1.0 
DRIMDA(1,2)=  0.0 
DRIMDA{1,3)=  0.5*SPINRZ(I)*DT 
DRIMDA(2,1)=  0.0 
DRIMDA(2,2)=  1.0 
DRIMDA{2,3)=  0.0 
DRIMDA(3,1)=  -0.5*SPINRZ(I) *DT 
DRIMDA(3,2)=  0.0 
DRIMDA(3,3)=  1.0 
CALL  INVERS(DRIMDA.3) 

ELSE 


LDWMDA(l,l)=(EXDOT(I)+ENSUM{I) ) 
LDWMDA(1,2)={EXYDOT(I)*0.5)  -  ZSPIN(I) 
LDWMDA{1,3)={EXZDOT{I)*0.5)  +  YSPIN(I) 
LDWMDA(2,1)=(EXYDOT{I)*0.5)  +  ZSPIN{I) 
LDWMDA(2,2)=(EYDOT(I)+ENSUM{I) ) 
LDWMDA(2,3)=(EYZDOT{I)*0.5)  -  XSPIN(I) 
LDWMDA(3,1)=(EXZDOT(I)*0.5)  -  YSPIN(I) 
LDWMDA(3,2)=(EYZDOT(I)*0.5)  +  XSPIN{I) 
LDWMDA(3,3)=(EZDOT(I)+ENSUM(I) ) 


DETMDA { 1 , 1 ) = { EXDOT ( I ) +ENSUM ( I ) ) *DT 
DETMDA ( 1 , 2 ) = ( EXYDOT ( I ) * 0 . 5 ) *DT 
DETMDA (1,3)= (EXZDOT (I ) *0.5) *DT 
DETMDA(2,1)= (EXYDOT (I) *0.5) *DT 
DETMDA(2,2)=(EYDOT(I)+ENSUM(I) ) *DT 
DETMDA(2,3)=(EYZDOT(I)*0.5)*DT 
DETMDA(3,1)=(EXZDOT(I)*0.5)*DT 
DETMDA(3,2)=(EYZDOT(I)*0.5) *DT 
DETMDA(3,3)=(EZDOT(I)+ENSUM(I) ) *DT 


DRMDA  (1,1)  = 
DRMDA  (1,2)  = 
DRMDA (1,3)= 
DRMDA (2,1)= 
DRMDA(2,2)= 
DRMDA(2,3)= 
DRMDA (3,1)= 
DRMDA(3,2)= 
DRMDA(3,3)= 


1.0 

-0.5*ZSPIN(I)*DT 
0.5*YSPIN(I)*DT 
0.5*ZSPIN(1) *DT 
1.0 

-0.5*XSPIN(I)*DT 

-0.5*YSPIN(I)*DT 

0.5*XSPIN(I)*DT 

1.0 


DRIMDA(1,1)=  1.0 
DRIMDA(1,2)=  0.5*ZSPIN(I)*DT 
DRIMDA(1,3)=  -0.5*YSPIN(I)*DT 
DRIMDA(2,1)=  -0.5*ZSPIN(I)*DT 
DRIMDA(2,2)=  1.0 
DRIMDA(2,3)=  0 . 5*XSPIN(I) *DT 
DRIMDA(3,1)=  0.5*YSPIN(I)*DT 
DRIMDA(3,2)=  -0 . 5*XSPIN(I) *DT 
DRIMDA(3,3)=  1.0 
CALL  INVERS(DRIMDA,3) 

END  IF 


C 
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RECOVER  CURRENT  DEFORMATION  GRADIENT 


FDGMDA (1,1) 
FDGMDA(1,2) 
FDGMDA (1,3) 
FDGMDA(2,1) 
FDGMDA(2,2) 
FDGMDA (2, 3) 
FDGMDA (3,1) 
FDGMDA(3,2) 
FDGMDA(3,3) 


=VHMDAV(1) 
=VHMDAV(2) 
=VHMDAV ( 3 ) 
=VHMDAV(4) 
=VHMDAV(5) 
=VHMDAV(6) 
=VHMDAV{7) 
=VHMDAV(8) 
=VHMDAV ( 9 ) 


ROTMDAd,!) 
ROTMDA (1,2) 
ROTMDAd,  3) 
R0TMDA(2,1) 
ROTMDA(2,2) 
ROTMDA(2,3) 
R0TMDA(3,1) 
ROTMDA(3,2) 
ROTMDA(3,3) 


=VHMDAV(20) 

=VHMDAV(21) 

=VHMDAV(22) 

=VHMDAV(23) 

=VHMDAV(24) 

=VHMDAV(25) 

=VHMDAV(26) 

=VHMDAV(27) 

=VHMDAV(28) 


CALL  REMDA ( DRMDA , DRIMDA , ROTMDA ) 


CONVERT  EPIC  (CAUCHY)  STRAIN  TO  GREEN- LAGRANGE  STRAIN 


CALC  DF  TO  SECOND  ORDER  TERM  LEVEL 
DO  IIN=1,3 
DO  JIN=1,3 
AUX=0.0 
DO  KIN=1,3 

AUX=AUX+LDWMDA ( IIN, KIN) *FDGMDA (KIN, JIN) 

END  DO 

DLF ( IIN, JIN) =AUX*DT 
END  DO 
END  DO 
DO  IIN=1,3 
DO  JIN=1,3 

AUX=DLF(IIN, JIN) 

DO  KIN=1,3 

AUX=AUX+DT2*LDWMDA(IIN, KIN) *DLF (KIN, JIN) 

END  DO 

DFdIN,  JIN)=AUX 
END  DO 
END  DO 

Calculate  Green-Lagrange  finite  strain  E  and  DE 
(only  the  upper  triangular  part) 

DO  IIN=1,3 
DO  JIN=1,3 
AUX1=0 . 0 
AUX2=0.0 
DO  KIN=1,3 

AUX1=AUX1+FDGMDA(KIN, IIN) *FDGMDA(KIN, JIN) 

AUX2=AUX2+DF (KIN, IIN) *FDGMDA(KIN, JIN) 

+FDGMDA(KIN, IIN) *DF (KIN, JIN) +DF (KIN, IIN) *DF (KIN, JIN) 


1 
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END  DO 

EIN { IIN, JIN) =AUX1 /2 . 0 
DEIN (IIN, JIN) =AUX2/2 . 0 
END  DO 

EIN (IIN, IIN) =EIN ( IIN, IIN) -0.50 
END  DO 


Extract  the  strains  and  strain  increments  into  col;ainn  matrices 
DO  IIN=1,3 

EJDC (IIN) =EIN ( IIN, IIN) +DEIN ( IIN, IIN) 

END  DO 

EJDC ( 4 ) =EIN (1,2) +DEIN (1,2) 

EJDC(6)=EIN(1,3)+DEIN(1,3) 

EJDC(5)=EIN(2,3)+DEIN(2,3) 

Update  the  deformation  gradient  F 
DO  JIN=1,3 
DO  IIN=1,3 

FDGMDA (IIN, JIN) =FDGMDA( IIN, JIN) +DF ( IIN, JIN) 

END  DO 
END  DO 


CONVERT  EJDC(I.NE.J)  TO  GAMMA 

DO  120  K=4,6,l 

EJDC (K)=EJDC(K) *2.0 
120  CONTINUE 

NONLINEAR,  INELASTIC  FRACTURE  MODEL 
EJDC  -  strain 

sigmda  -  stress  returned  by  the  model 
There  are  five  history  quantities. 

DO  IHyPO=10,15,l 

HISJDC(IHYPO-9)  =  VHMDAV(IHYPO) 

END  DO 

CALCULATE  VOLUMETRIC  STRAIN 

DVOL  =  VOLUMETRIC  STRAIN  CALCULATED  BY  EPIC;  TENSION  + 

EPSKK  =  -DVOL (I) 

EPSKK  =  -EJDC(l)  -  EJDC(2)  -  EJDC(3) 

ASSIGN  HISTORY  VARIABLES 

EPSKKM  =  HISJDC(l) 

ATYCP  =  HISJDC(2) 

AGAMCP  =  HISJDC(3) 

GMOCTM  =  HISJDC(4) 

FLGJDC  =  HISJDC(5) 

INITIALIZE  SOME  VARIABLES 

EPSKKO  =0.0 
AKULRL  =  AK 
COHSN  =0.0 


oono  oo  o  n  oooooo  noo  on  non 


194 


FRIC 

= 

0.0 

DEV 

_ 

0.0 

SIGM 

z= 

0.0 

SIGMG 

= 

0.0 

FACTRF 

= 

1.325 

ATYC 

= 

0.0 

AGAMC 

= 

0.005 

EBARP 

= 

0.0 

EBARS 

=: 

0.0 

GMOD 

= 

0.0 

GULRL 

=: 

400000.0  +  AMUC 

BLKJDC 

AK 

CALCULATE  OCTAHEDRAL  SHEAR  STRAIN 

E1E2  =  (EJDCd)  -  EJDC{2)  )**2.0 
E1E3  =  (EJDCd)  -  EJDC(3))**2.0 
E2E3  =  (EJDC(2)  -  EJDC (3 ) ) **2 . 0 
EPS12  =  EJDC (4)  *  EJDC (4) 

EPS13  =  EJDC{6)  *  EJDC(6) 

EPS23  =  EJDC (5)  *  EJDC (5) 

****  use  1.5  instead  of  6.0  since  eps  =  gamma  /  2  and  current 
****  values  of  EJDC (4-6)  are  gamma  and  this  eqn  assumes  eps 
GAMOCT  =  TW03RD  *  SQRT(E1E2  +  E1E3  +  E2E3 
1  +  1.5  *  (EPS12  +  EPS13  +  EPS23)) 

HYDROSTATIC  PART 

EPSKKV  =  AMAXl (EPSKK,EPSKKM) 

HYDROSTATIC  COMPRESSION  PART 


LOADING  IN  COMPRESSION 

IF  ( EPSKKV. LT.AEPSC)  THEN 
SIGM  =  AK  *  EPSKKV 
BLKJDC  =  AK 

ELSEIF  { EPSKKV. GT.AEPSL)  THEN 

SIGM  =  SIGML  +  (EPSKKV  -  AEPSL)  *  AKL 
BLKJDC  =  AKL 

ELSE 

EPSKKP  =  EPSKKV  -  AEPSC 

SIGM  =  AKEPSC  +  AAH  *  EPSKKP  +  ABH  *  EPSKKP*EPSKKP 
&  +  ACH  *  EPSKKP* *3.0 

BLKJDC  =  AK  +  AKLAK*EPSKKP/AEPSLC 


ENDIF 


UNLOAD /RELOAD  IN  COMPRESSION 


IF  (EPSKKM.GT.EPSKK)  THEN 
IF  (EPSKKM.lt. AEPSC)  THEN 
AKULRL  =  AK 

ELSEIF  (EPSKKM. GT.AEPSL)  THEN 
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AKULRL  =  AKL 
ELSE 

AKULRL  =  AK  +  AKLAK* (EPSKKM  -  AEPSC ) / AEPSLC 
END  IF 

EPSKKO  =  EPSKKM  -  (SIGM/ AKULRL) 

EPSKKT  =  EPSKKO  -  (ASIGS /AKULRL) 

SIGM  =  SIGM  -  AKULRL  *  (EPSKKM  -  EPSKK) 

BLKJDC  =  AKULRL 
IF  (EPSKK.lt. EPSKKT)  THEN 
EPSKKT  =  EPSKKT  -  EPSKK 
ELSE 

EPSKKT  =0.0 
END  IF 

IF  ( SIGM. LT. -ASIGS)  SIGM  =  -ASIGS 
EBARP  =  1.0  *  EPSKKT  /  AEPSS 
END  IF 
C 

C  LOADING  IN  TENSION 

C 

ca  IF  (EPSKK.lt. 0.0  .OR.  SIGM. LT. 0.0)  THEN 

C 

C  EPSKKO  IS  STRAIN  AT  PRESSURE  =  ZERO 

C 

ca  EPSKKT  =  ABS (EPSKK  -  EPSKKO) 

ca  RNUM  =  AKULRL  *  (ASIGS  /  ASIGF  -  1.0) 

ca  RDEN  =  ASIGS/AEPSS  *  (AEPSF/AEPSS  -  1.0)  **  2.0 

ca  R  =  RNUM  /  RDEN  -  AEPSS /AEPSF 

ca  A  =  R  +  AKULRL* AEPSS /ASIGS  -  2.0 

ca  B=2.0*R-1.0 

ca  X  =  EPSKKT  /  AEPSS 

ca  SIGM  =  -  AKULRL  *  EPSKKT  /  (1.0  +  A*X  -  B*X*X  +  R*X*X*X) 

C 

C  NO  SOFTENING  WITH  FOLLOWING  LINE 

ca  IF  (EPSKKT. GT. AEPSS)  SIGM  =  -ASIGS 

ca  IF  ( SIGM. LT. -ASIGS)  SIGM  =  -ASIGS 

ca  BLKJDC  =  AKULRL 

c 

c  Check  for  failure  in  tension 

ca  EBARP  =  EPSKKT  /  AEPSS 

C 

ca  END  IF 

C 

IF  (GAMOCT.GT.GMOCTM)  THEN 
SIGMG  =  SIGM 
ELSE 

SIGMG  =  HISJDC(6) 

END  IF 

C  SIGMG  =  SIGM 

C 

C  BEGIN  DEVIATORIC  PART 

C 

GAM  =  AMAXl (GAMOCT,GMOCTM) 

C 

C  LOADING 

C 

C  FRICTION  SUBPART 

c 

c  Calculate  ultimate  surface 

C 


o  n 
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c 

IF  (SIGMG.GT.0.0)  THEN 

FALULT  =  AA  -  AC  *  EXP{AB  *  SIGMG) 

ELSE 

FALULT  =  AAAC  *  (1.0  +  SIGMG/ASIGS) 

END  IF 
C 

IF  (TIME.LE.DT)  ATYCP  =  AAAC  *  FACTRF 
c 
c 

c  Calculate  cohesion  yield  stress;  to  be  subtracted  from 

c  FALULT  to  give  friction  yield  stress 

c 

c  Increase  yield  stresses  by  FACTRF  to  account  for  hyperbolic 

c 

IF  (SIGMG.GT.PFRIC)  THEN 
FACTRF  =1.1 
ATYC  =0.0 
ELSE 

IF  (SIGMG.lt. 0.0)  THEN 
ATYC  =  FALULT 
FACTRF  =  1.325 

ELSEIF  (SIGMG.GT.PCRUSH)  THEN 

ATYC  =  AAAC* (1.0- (SIGMG  -  PCRUSH ) / ( PFRIC  -  PCRUSH) ) 

FACTRF  =  1.1  +  0.225  *  (1.0  -  (SIGMG  -  PCRUSH) 

1  /(PFRIC  -  PCRUSH)) 

ELSE 

ATYC  =  AAAC 
FACTRF  =  1.325 
END  IF 
END  IF 

ATYC  =  ATYC  *  FACTRF 
ATYCP  =  AMINl (ATYCP, ATYC) 

FALJDC  =  FALULT  *  FACTRF 
c 

c  Calculate  friction  subpart 

C  AND 

c  Calculate  agamc  at  fail(fric) 

c 
c 

IF (SIGMG.GT.0.0)  THEN 
C 

c  Can  have  friction  with  sigm  <  0.0 

c 

GMOD  =  400000.0  +  ABO  *  SIGMG**ABl 

AGAMC  =  FACTRF* (FALJDC  -  ATYCP) * (FALJDC/ FACTRF  -  ATYCP)/ 

1  AMAX1(ABS( SMALL), (GMOD* FALJDC* (FACTRF  -  1.0))) 

ELSE 

GMOD  =  400000.0 
AGAMC  =  0.005 
END  IF 

Limit  AGAMC  to  be  greater  than  0.005  (0.5%) 

AGAMC  =  AMAXl (AGAMC, AGAMC P, 0.005) 
c 

FRICDEN  =  FALJDC  -  ATYCP  +  GAM  *  GMOD 
IF  (FALJDC-ATYCP.LT.ABS (SMALL) )  THEN 
FRIC  =0.0 
ELSE 
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FRIC  =  (FALJDC  -  ATYCP)  *  GMOD  /  FRICDEN 
END  IF 
c 

C  COHESION  SUBPART 

C 

C 

IF  (ATYCP.GT.ABS (SMALL) )  THEN 
REDUC  =1.0 
c 

c  Reduction  in  cohesion  subpart 

c 

IF  (GAM.GT.AGAMC)  THEN 

REDUC  =  1.0  -  (GAM  -  AGAMC) / (AGAMC* (ANGAMC  -  1.0)) 

END  IF 

REDUC  =  AMAXl (REDUC, 0.0) 

C 

C  NO  SOFTENING  WITH  FOLLOWING  LINE 

C  REDUC  =1.0 

C 

c  Calculate  cohesion  subpart  prior  to  reduction 

c 

COHSN  =  ATYCP  *  AGAMC  *  AMUC  *  REDUC  / 

1  (ATYCP*AGAMC  +  (AMUC*AGAMC  -  ATYCP) *GAM) 

ELSE 

COHSN  =0.0 
END  IF 
C 

c  Add  friction  and  cohesion  siabparts 

c 

DEV  =  FRIC  +  COHSN 
IF  (SIGM.LE.-ASIGS)  DEV  =  0.0 
C  GULRL  =  DEV 

IF  (SIGMG.LE.0.0)  THEN 
GULRL  =  400000.0 
ELSE 

GULRL  =  400000.0  +  AB0*SIGMG**AB1 
END  IF 

GULRL  =  GULRL  +  AMUC 
C 

C  UNLOAD /RELOAD 

C 

ca  REMOVE  ca  COMMENT  CHARACTERS  TO  RETURN  TO  USE  NONLINEAR 
ca  UNLOADING /RELOADING 
C 

IF  (GAMOCT.LT.GMOCTM)  THEN 
TOCTM  =  DEV  *  GMOCTM 

ca  IF  (TOCTM.lt. 0.001*FALJDC)  THEN 

ca  TAUOCT  =  TOCTM 

ca  ELSE 

c  Initial  uloading  modulus  based  on  pressure  at  start  of  tinloading 

IF  (SIGMG.LE.0.0)  THEN 
GULRL  =  400000.0 
ELSE 

GULRL  =  400000.0  +  AB0*SIGMG**AB1 
END  IF 

GULRL  =  GULRL  +  AMUC 
GULRL  =  GULRL  *  1.20 

ca  GMOCTO  =  GMOCTM  -  TOCTM  /  GULRL 

ca  IF  (GAMOCT.GE. GMOCTO)  THEN 


198 


ca 

ca 


c 

c 

c 

c 

c 

c 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

c 

c 

C 

c 

c 

ca 

c 

c 

c 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

ca 

c 

c 

ca 

C 

ca 

ca 

ca 

ca 

ca 


TAUOCT  =  ABS{TOCTM  -  (GMOCTM  -  GAMOCT) *GULRL) 

ELSE 

GAM  =  GMOCTO  -  GAMOCT 
IF  (SIGM.GT.0.0)  THEN 

FALULT  =  AA  -  AC  *  EXP(AB  *  SIGM) 

ELSE 

FALULT  =  AAAC  *  (1.0  +  SIGM/ASIGS) 

END  IF 

Calculate  cohesion  yield  stress;  to  be  subtracted  from 
FALULT  to  give  friction  yield  stress 

Increase  yield  stresses  by  FACTRF  to  accoiint  for  hyperbolic 

IF  (SIGM.GT.PFRIC)  THEN 
FACTRF  =1.1 
ATYCUR  =0.0 
ELSE 

IF  (SIGM. LT. 0.0)  THEN 
ATYCUR  =  FALULT 
FACTRF  =  1.325 

ELSEIF  (SIGM.GT.PCRUSH)  THEN 

ATYCUR  =  AAAC* (1.0- (SIGM-PCRUSH) / (PFRIC-PCRUSH) ) 
FACTRF  =  1.1  +  0.225  *  (1.0  -  (SIGM  -  PCRUSH) 

1  /(PFRIC  -  PCRUSH)) 

ELSE 

ATYCUR  =  AAAC 
FACTRF  =  1.325 
END  IF 
END  IF 

ATYCUR  =  ATYCUR  *  FACTRF 
ATYCPUR  =  AMINKATYCP,  ATYCUR) 

FALJDC  =  FALULT  *  FACTRF 

Calculate  friction  subpart 
AND 

Calculate  agamcur  at  fail(fric) 

IF ( SIGM. LE. 0.0)  THEN 

Can  have  friction  with  sigm  <  0.0 

GMOD  =  400000.0 
AGAMCUR  =  0.005 
ELiSS 

GMOD  =  400000.0  +  ABO  *  SIGM**AB1 
AGAMCUR  =  FACTRF* (FALJDC  -  ATYCPUR) 

1  *( FALJDC /FACTRF  -  ATYCPUR)/ 

2  AMAXl (ABS (SMALL) , (GMOD*FALJDC* (FACTRF  -  1.0))) 

END  IF 

Limit  AGAMCUR  to  be  greater  than  0.005  (0.5%) 

AGAMCUR  =  AMAXl (AGAMCUR, AGAMCP, 0.005) 

FRICDEN  =  FALJDC  -  ATYCPUR  +  GAM  *  GMOD 
IF  ( FALJDC- ATYCPUR. LE. 0.0)  THEN 
FRIC  =0.0 
ELSE 

FRIC  =  (FALJDC  -  ATYCPUR)  *  GMOD  /  FRICDEN 
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ca 

c 

C 

c 

ca 

c 

C 

C 

ca 

c 

c 

c 

ca 

ca 

ca 

ca 

ca 

C 

c 

c 

ca 

ca 

ca 

ca 

ca 

ca 

ca 


c 

C 

C 

C 

c 

c 

c 


c 


END  IF 

COHESION  SUBPART 
IF  (ATYCPUR.GT.0.0)  THEN 

NO  SOFTENING  WITH  FOLLOWING  LINE  ALLOWED  IN  UNLOAD/RELOAD 
REDUC  =1.0 

Calculate  cohesion  subpart  prior  to  reduction 

COHSN  =  ATYCPUR  *  AGAMCUR  *  AMUC  *  REDUC 
1  /  (ATYCPUR*AGAMCUR  +  (AMUC*AGAMCUR  -  ATYCPUR) * GAM) 

ELSE 

COHSN  =0.0 
END  IF 

Add  friction  and  cohesion  subparts 

TAUOCT  =  AMIN1{{FRIC  +  COHSN)  *  GAM,FALULT) 

END  IF 
END  IF 

DEVU  =  -TAUOCT/GAMOCT 
DEV  =  DEVU 
DEV  =  ABS(DEV) 

Following  line  added  with  ca  changes 
IF  (TAUOCT. GT. 0.0)  THEN 

TAUOCT  =  AMINl (TAUOCT, FALULT) 

ELSE 

TAUOCT  =  AMAXl (TAUOCT, -FALULT) 

END  IF 

DEV  =  TAUOCT/GAMOCT 
IF  (SIGM.LE.-ASIGS)  DEV  =  0.0 

END  IF 

STRESS  CALCUALTION 

EBARS  =  GAMOCT  /  AGAMC 
If  EBARP  or  EBARS  GT  EFAILV  then 
set  EBARS  and  EBARP  to  be  big 
allow  pressure  only,  i.e.,  DEV=0.0 
IF  ( AMAXl ( EBARP , EBARS ) . GT . EFAILV )  THEN 
EBARS  =  EBARS*100.0 
EBARP  =  EBARP*100.0 
DEV  =0.0 
ENDIF 

IF  (HMIN(I) .LT.0.05)  THEN 
EBARS  =  EFAILV*100.0 
EBARP  =  EFAILV*100.0 
ENDIF 

EPSKKD  =  EJDC(l)  +  EJDC(2)  +  EJDC(3) 

DEVEPS  =  DEV  *  EPSKKD  *  THIRD 

SIGMDA(l)  =  -SIGM  +  2 . 0* (DEV  *  EJDC(l)  -  DEVEPS) 

SIGMDA(2)  =  -SIGM  +  2.0* (DEV  *  EJDC(2)  -  DEVEPS) 

SIGMDAO)  =  -SIGM  +  2.0*(DEV  *  EJDC(3)  -  DEVEPS) 

VHMDA(J,19)  =  (SIGMDA(l)  +  SIGMDA(2)  +  SIGMDA(3))  *  THIRD 

SIGMDA(4)  =  DEV  *  EJDC(4) 


oo  oooooooooo  o  o  oooo 
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SIGMDA(5)  =  DEV  *  EJDC(5) 

SIGMDA{6)  =  DEV  *  EJDC(6) 

HISJDC(7)  =  SIGMG 

HISJDC(8)  =  DEV 


UPDATE  HISTORY  PARAMETERS 

HISJDCd)  =  EPSKKV 

HISJDC(2)  =  ATYCP 

HISJDCO)  =  AGAMC 

HISJDC{4)  =  AMAXl(GAMOCT,GMOCTM) 

HISJDC(5)  =  AMAX1{EBARP,EBARS) 

EBAR(I)  =  HISJDC(5) 

HISJDC(6)  =  SIGMG 

VHMDA(J,29)  =  EJDC(l) 

VHMDA(J,30)  =  EJDC{2) 

VHMDA(J,31)  =  EJDC{3) 

****  convert  gamma  back  to  eps  **** 
VHMDA(J,32)  =  EJDC(4)  /  2.0 
VHMDA(J,33)  =  EJDC(5)  /  2.0 
VHMDA(J,34)  =  EJDC(6)  /  2.0 
VHMDA{J,35)  =  SIGM 
VHMDA(J,36)  =  SIGMDA(l) 
VHMDA(J,37)  =  SIGMDA{2) 
VHMDA(J,38)  =  SIGMDA(3) 
VHMDA(J,39)  =  SIGMDA{4) 
VHMDA{J,40)  =  SIGMDA(5) 
VHMDA(J,41)  =  SIGMDA(6) 


Rotate  stresses  back  to  current  orientation 
CALL  RSMDA{ROTMDA,SIGMDA) 

CALCULATE  "4G/3"  FOR  SOUND  SPEED  CALCULATION 
CGVJDC  =  4.0  *  GULRL  *  THIRD 
SOUND  SPEED 

CMODJDC  =  1.0*(BLKJDC  +  CGVJDC) 

BULK { I )  =  BLK JDC 

SS2 ( I ) =ABS (CMODJDC/ (DENV/ ( 1 . +DVOL ( I ) ) ) ) 

SS2(I)  =  AMAX1(SS2 (I) ,1.0) 

!  TOTAL  STRESSES  ! 

SX(I)  =  SIGMDAd) 

SYd)  =  SIGMDA{2) 

SZd)  =  SIGMDA(3) 

SXY(I)  =  SIGMDA(4) 

SYZd)  =  SIGMDA(5) 

SXZ(I)  =  SIGMDA(6) 

IF  (IGEOM.lt. 7)  THEN 
SYZ(I)  =0.0 
SXY(I)  =  0.0 
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c 

c 

c 

c 


C 

C 

C 


C 


C 


130 

131 
c 

c 

c 


c 


C 

C 

100 

C 

C 

C 

C 


C 

C 


ENDIF 

SBMDA(I)  =  (SX(I)  +  SY(I)  +  SZ(I))*THIRD 
CONVERT  TO  "DEVIATORIC"  STRESS 


SX(I)  =  SX{I)  -  SBMDA(I) 

SY(I)  =  SY{I)  -  SBMDA(I) 

SZ(I)  =  SZ(I)  -  SBMDA(I) 

STORE  UPDATED  DEFORMATION  GRADIENT 


VHMDAV(l) 
VHMDAV{2) 
VHMDAV(3) 
VHMDAV(4) 
VHMDAV{5) 
VHMDAV(6) 
VHMDAV{7) 
VHMDAV{8) 
VHMDAV ( 9 ) 


FDGMDA(1,1) 

FDGMDA(1,2) 

FDGMDA{1,3) 

FDGMDA{2,1) 

FDGMDA{2,2) 

FDGMDA(2,3) 

FDGMDA(3,1) 

FDGMDA(3,2) 

FDGMDA(3,3) 


VHMDAV (20) 
VHMDAV(21) 
VHMDAV (22) 
VHMDAV (23) 
VHMDAV (24) 
VHMDAV(25) 
VHMDAV(26) 
VHMDAV(27) 
VHMDAV(28) 


ROTMDAd,!) 
ROTMDA(l,2) 
ROTMDAd,  3) 
ROTMDA(2,l) 
ROTMDA(2,2) 
ROTMDA(2,3) 
ROTMDA(3,l) 
ROTMDA(3,2) 
ROTMDA(3,3) 


DO  130  K=l,9,l 

VHMDA ( J , K )= VHMDAV ( K ) 

CONTINUE 

DO  131  K=20,28,l 

VHMDA ( J , K )= VHMDAV ( K ) 

CONTINUE 

*************** 

changed  from  10,14,1  to  10,13,1  for  corant  check  for  failure 

DO  IHYPO=10,15,1 

VHMDA ( J , IHYPO ) =HIS JDC ( IHYPO- 9 ) 

END  DO 

★*★*★**★★***★** 

VHMDA (J, 16)  =  GULRL  *  DT 
VHMDA(J,17)  =  HISJDC(8) 

VHMDA(J,18)  =  HISJDC(5) 

****★★*★★***★*★ 


CONTINUE 


TEMPORARY  APPROXIMATION 

IF ( IGEOM . LE . 7 ) THEN 
1  =  0 

DIR$  VECTOR 
VD$  VECTOR 
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DO  150, 

I  =  I 
EDEV(I) 


J=Ll,LlSr 
+  1 


1 

2 

3 

4 

5 

6 


(SX(I)+SX1(I) )*EXDOT(I) 
(SZ{I)+S21(I) )*EZDOT(I) 
(SY{I)+SY1(I) )*EYDOT(I) 
(SXZ(I)+SXZ1(I) )*EXZDOT(I) 
(SXY(I)+SXY1(I) )*EXYDOT{I) 
{SYZ{I)+SYZ1(I) )*EYZDOT(I) ) 


EXDOT (I) 

EZDOT ( I ) 

EYDOTd) 

EXZDOTd) 

EXYDOTd) 

EYZDOTd) 


(l.O+DVBARd)  )  *DT2 


EXDOT  d) 
EZDOT (I ) 
EYDOT (I ) 
EXZDOTd) 
EXYDOTd) 
EYZDOTd) 


1 

2 

3 

4 

5 


1 

2 

3 

4 

5 

6 


IFdCHECKd)  .GE.0)THEN 
EPDOTd)  =  SQRT(C2D9 


(SXd)  -  SXld) 

(szd)  -  szid) 

(SYd)  -  SYld) 

■  (sxzd)-sxzid) 

•  (SXYd)-SXYld) 

■  {SYZd)-SYZld) 


-  DSXd))*TGDTI 

-  DSZd))*TGDTI 

) *TGDTI 

-DSXZd)  )  *GDTI 
) *GDTI 
) *GDTI 


'(  (EXDOTd) -EZDOTd)  )  **2 
+  (EXDOTd) -EYDOTd)  )  **2 
+  (EZDOTd) -EYDOTd)  )**2 
+  1.5* (  EXZDOT(I)*EXZDOT(I) 

+  EXYDOT(I)*EXYDOT(I) 

+  EYZDOT ( I ) *EYZDOT ( I ) ) ) ) 
EBAR(I)  =  EBAR(I)  +  EPDOT(I)*DT 
EP(I)  =  EP(I)  +  (  (SX(I)+SX1(I) )*EXDOT(I) 

(SZ(I)+SZ1(I) )*EZDOT(I) 
(SY(I)+SY1(I) )*EYDOT(I) 
(SXZ(I)+SXZld)  )*EXZDOTd) 
(SXY(I)+SXY1 (I) ) *EXYDOT(I) 
(SYZ(I)+SYZ1(I) )*EYZDOT(I) ) 
(l.O+DVBARd)  )*DT2 


+ 

+ 

+ 

+ 


ENDIF 

150  CONTINUE 

enddo 

IF ( IGEOM . EQ . 4 ) THEN 
INCOMP  =  1 
ENDIF 

IF (IGEOM. GE. 5 
INCOMP  =  2 
EIODIF 
ELSE 

here  3-D 
1  =  0 

DO  170,  J=L1,LN 
1  =  1  +  1 
EDEV(I) 


AND.  NODE3 . EQ . 0 ) THEN 


1 

2 

3 

4 

5 

6 


{ 

+ 

+ 

+ 

+ 

+ 


(SX(I)+SX1(I) ) *EXDOT(I) 

(SY(I)+SY1(I) )*EYDOT(I) 

(SZ(I)+SZ1(I) )*EZDOT(I) 

(SXY(I)+SXY1(I) )*EXYDOT(I) 
(SXZ(I)+SXZ1(I) )*EXZDOT(I) 
(SYZ(I)+SYZ1(I) )*EYZDOT(I) ) 

*  (l.O+DVBARd)  )*DT2 

EXDOT(I)  =  EXDOT (I)  -  (SX(I) -SXl (I) -DSX(I) ) *TGDTI 
EYDOT(I)  =  EYDOTd)  -  (SY(I) -SYl  (I) -DSY(I)  )  *TGDTI 
EZDOT(I)  =  EZDOT (I)  -  (SZ ( I) -SZl ( I ) -DSZ (I ) ) *TGDTI 
EXYDOT(I)  =  EXYDOTd )  -  (SXY(I) -SXYl  (I) -DSXY(I)  )  *GDTI 
EXZDOT(I)  =  EXZDOTd )  -  (SXZ  (I) -SXZl  (I) -DSXZ  (I)  )  *GDTI 
EYZDOT (I)  =  EYZDOT (I)  -  (SYZ (I) -SYZl (I) -DSYZ (I) ) *GDTI 
IFdCHECKd)  .GE.0)THEN 

EPDOT(I)  =  SQRT(C2D9*(  (EXDOT(I) -EYDOT(I) ) **2 
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1 

2 

3 

4 

5 
c 

1 

2 

3 

4 

5 

6 

170 

C 


+  (EXDOT(I)-EZDOT{I) )**2 
+  (EYDOTd) -EZDOT(I)  )  **2 
+  1.5*  (  EXYDOT(I)*EXyDOT(I) 

+  EXZDOT(I) *EXZDOT(I) 

+  EYZDOT(I) *EYZDOT(I) ) ) ) 
EBAR(I)  =  EBARd)  +  EPDOTd)*DT 
EPd)  =  EPd)  +  (  (SXd)+SXld)  )*EXDOTd) 

+  (SYd)+SYld)  )*EYDOTd) 

+  (SZd)+SZld)  )*EZDOTd) 

+  (SXYd)+SXYl  (I)  )  *EXYDOTd) 

+  (SXZd)+SXZl  d)  )  *EXZDOTd) 

+  (SYZd)+SYZld)  )*EYZDOTd)  ) 

*  d.O+DVBARd)  )  *DT2 

ENDIF 

CONTINUE 

enddo 

IF(NODE3.EQ.O)THEN 
INCOMP  =  3 

ELSEIF (NODE4 . EQ . 0 ) THEN 
INCOMP  =  4 
ENDIF 


ENDIF 


c 

c  Convert  to  total  stresses 


C  DIR$  NOVECTOR 

1  =  0 

DO  200,  J=L1,LN 
1  =  1  +  1 

IFdCHECKd)  .LT,0)GOTO  200 
sxd)  =  SXd)  +  SBMDAd) 

SYd)  =  SYd)  +  SBMDAd) 

SZd)  =  SZd)  +  SBMDAd) 

SBARd)  =  SBMDAd) 

200  CONTINUE 


C 


RETURN 


END 

*$ 

SUBROUTINE  INVERS(D,NE) 

Q  ******************************************************************** 

c  *  * 

C  *  INVERS  -  TO  COMPUTE  THE  INVERSE  OF  A  MATRIX  USING  * 

C  *  GAUSS  ELIMINATION  * 

C  *  * 

Q  ******************************************************************** 

REAL  CI(3,6) ,D(3,3) 

REAL  TMIJ,SUM 
INTEGER  NE,I,J,NN,K 
CDIR$  NOVECTOR 
C 

DO  3  I=1,NE,1 
DO  5  J=1,NE,1 
CId,J)=Dd,J) 

Dd,  J)=0.0 
Cld,  J+NE)=0.0 
5  CONTINUE 

CId,NE+I)=1.0 
3  CONTINUE 
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c 

Q  ***  START  THE  GAUSS  ELIMINATION  PROCEDURE  *** 

C 

NN=2*NE 

DO  10  J=1,NE-1,1 
DO  40  I=J+1,NE,1 

TMIJ=CI ( I , J) /Cl ( J, J) 

DO  50  K=J+1,NN,1 

CI(I,K)=CI{I,K)-TMIJ*CI(J,K) 

50  CONTINUE 

40  CONTINUE 

10  CONTINUE 
C 

C  ***  UTILIZE  BACK-SUBSTITUTION  TO  OBTAIN  THE  INVERSE  *** 

C 

DO  55  K=1,NE,1 
DO  60  I=NE,1,-1 
SUM=0 . 0 

DO  70  J=I+1,NE,1 

SUM=SUM+CI ( I , J ) *D ( J, K) 

70  CONTINUE 

D{I,K)={CI(I,NE+K)-SUM)/CI{I,I) 

60  CONTINUE 

55  CONTINUE 
C 
C 

999  RETURN 
END 

*$ 

^  ********************************************************************** 
C  ***  SUBROUTINE  REMDA  -  TO  ROTATE  STRAIN  INCREMENT  TO  INITIAL  CONFIG.  * 

Q  ********************************************************************** 

SUBROUTINE  REMDA { DRMDA , DRIMDA , ROTMDA ) 

REAL  DRMDA(3,3) ,DRIMDA(3,3) ,ROTMDA(3,3) 

REAL  DELTAR{3,3) ,T1(3,3),T2{3,3),C(3,3) 

INTEGER  I,J,K 
CDIR$  NOVECTOR 
C 
C 

C  DELTAR  =  DRIMDA*DRMDA 

C 

DO  1=1,3 
DO  J=l,3 

DELTAR ( I, J) =0.0 
DO  K=l,3 

DELTAR ( I , J) =DELTAR { I , J) +DRIMDA ( I , K) *DRMDA (K, J) 

END  DO 
END  DO 
END  DO 
C 

C  ROTMDA_T+l  =  DELTAR *ROTMDA_T 

C 

DO  1=1,3 
DO  J=l,3 

T1{I, J)=0.0 
DO  K=l,3 

T1 (I, J)=T1{I, J)+DELTAR{I,K) *ROTMDA(K, J) 

END  DO 
END  DO 
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END  DO 
C 

C  COMPUTE  ROTMDA_TRANSPOSE*DETMDA 

C 

DO  1=1,3 
DO  J=l,3 

ROTMDA(I,J)  =  T1(I,J) 

END  DO 
END  DO 
999  RETURN 
END 

*$ 

Q  ********************************************************************** 

C  ***  SUBROUTINE  RSMDA  -  TO  ROTATE  STRESSES  TO  CURRENT  CONFIGURATION  * 

Q  ********************************************************************** 

SUBROUTINE  RSMDA ( ROTMDA, SIGMDA) 

REAL  ROTMDA(3,3) ,T1(3,3),T2(3,3) ,SIGMDA(6) 

INTEGER  I,J,K 
CDIR$  NOVECTOR 
C 

C  insert  the  stresses  into  square  matrices 
C 

T1(1,1)=SIGMDA{1) 

T1(1,2)=SIGMDA{4) 

T1(1,3)=SIGMDA{6) 

T1(2,1)=SIGMDA{4) 

T1(2,2)=SIGMDA(2) 

T1(2,3)=SIGMDA{5) 

T1{3,1)=SIGMDA{6) 

T1(3,2)=SIGMDA(5) 

T1(3,3)=SIGMDA{3) 

C 

C  COMPUTE  ROTMDA* STRESS 

C 

DO  1=1,3 
DO  J=l,3 

T2(I,J)=0.0 
DO  K=l,3 

T2 (I, J)=T2 (I, J)+ROTMDA(I,K)*Tl(K, J) 

END  DO 
END  DO 
END  DO 
C 

C  COMPUTE  { ROTMDA* STRESS ) *ROTMDA_TRANSPOSE 

C 

DO  1=1,3 
DO  J=l,3 

T1{I,J)=0.0 
DO  K=l,3 

T1{I,  J)=T1(I,  J)+T2(I,K)*ROTMDA{J,K) 

END  DO 
END  DO 
END  DO 
C 

C  INSERT  THE  STRESSES  INTO  A  VECTOR 

C 

SIGMDA{1)=T1{1,1) 

SIGMDA{2)=T1(2,2) 

SIGMDA(3)=T1(3,3) 


SIGMDA(4)=T1{1,2) 
SIGMDA(5)=T1 (2,3) 
SIGMDA{6)=T1(1,3) 


999  RETURN 
END 
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